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1.0  StXMARY 


An  accurate  and  rapid  method  is  presented  for  solution  of  the  general 
equations  of  compressible,  steady,  laminar-boundary-layer  flow.  The  method 
allows  arbitrary  conditions  on  all  of  the  following:  pressure  gradient,  sur¬ 
face  temperature  and  its  gradient,  heat  transfer,  mass  transfer,  and  fluid 
properties.  Also,  the  method  can  calculate  the  second-order  effect  of  trans¬ 
verse  curvature .  The  only  restrictions  on  the  method  are  that  the  body  be 
axially  symmetric  or  two-dimensional  and  that  no  dissociation  of  the  fluid 
occurs. 

The  equations  that  are  solved  are  developed  from  the  Navier-Stokes  and 
energy  equations  by  an  order-of-magnitude  analysis.  They  differ  from  the  con¬ 
ventional  boundary-layer  equations  of  Prandtl  only  in  that  the  second-order 
terms  that  include  transverse  curvature  are  retained. 

The  method  of  solution  consists  of  replacing  the  partial  derivatives 
with  respect  to  the  flow  direction  by  finite  differences,  while  retaining  the 
derivatives  in  a  direction  normal  to  the  boundary,  so  that  the  partial  differen¬ 
tial  equations  become  approximated  by  ordinary  differential  equations.  Reasons 
for  choosing  this  method  rather  than  the  more  conventional  finite-difference 
methods  are  discussed. 

Arbitrary  fluid  properties  may  be  vised  in  the  method  of  solution,  that  is, 
they  are  inputs  in  the  computer  program  in  the  form  of  formulas  or  tables  as 
functions  of  local  enthalpy  and  pressure.  Results  obtained  with  the  method 
using  exact  fluid  properties  for  air  are  compared  with  those  using  simpler 
fluid-property  laws.  These  simpler  laws,  which  often  have  been  used  in  the 
past  in  boundary-layer  investigations,  ere  shown  often  to  give  poor  predic¬ 
tions  of  heat  transfer  and  skin  friction  at  high  speeds. 

The  method  has  been  programmed  on  the  IBM  7090  computer,  and  solutions 
for  a  wide  variety  of  flows  are  presented.  Comparisons  are  made  with  other 
exact  and  approximate  methods  of  solutions.  The  flows  include  cases  of  heat 
transfer,  mass  transfer,  and  discontinuities  in  the  boundary  conditions  over 
a  large  range  in  Mach  number  (Mach  0.0  to  10.0 ).  Some  comparison  with  ex- 
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perimental  measurement  is  also  made.  Also  a  study  of  the  effect  of  trans¬ 
verse  curvature  on  the  flow  over  cones  is  presented.  The  large  number  of 
calculations  and  comparisons  establish  that  the  method  is  rapid,  highly 
accurate,  and  powerful.  It  appears  capable  of  solving  any  flow  problem  for 
which  the  boundary-layer  equations  themselves  remain  valid. 
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4.0  PRINCIPAL  NOTATION 


local  velocity  of  sound 
constants  defined  vhere  applicable 


constant  coefficients  defined  where  applicable 

chord  length 

local  skin-friction  coefficient,  eq. (6.114) 

skin-friction  parameter,  eq.(6.1l6) 
specific  heat  at  constant  pressure 

£J±_ 

Pe“e 

nondimens ional  stream  function,  defined  by  eq.(6.7) 
body  force 

nondimens  ional  total  enthalpy  ratio  =>  H/H^ 
local  enthalpy 

reference  enthalpy  used  in  fluid  property  relation  from 
Ref.  7,  =  2.119  x  10®  ft2/sec2 

1  2 

total  enthalpy  =  h  +  -g  u 

count  of  successive  tries  in  procedure  for  solving  momentum 
equation  (Appendix  B) 
unit  vector,  eq.(A-5) 

thermal  conductivity 

bound  in  computer  program  on  values  of  q>* 

count  of  successive  solutions  of  the  energy  equation  and 
fluid  properties  in  the  procedure  of  solution 

the  exponent  in  the  free  stream  velocity  variation,  U  =  c^ 


M 

n 

N 

P 

P 

Pr 

q. 

Q, 


R 

Re 

St 

T 

u 

v 

V 

X 

X 

X 

y 

y 

y 


Mach  number 

count  of  the  number  of  steps  In  the  x-direction 

^+R 

pressure 

du 

x  e 

u  dx 
e 

Prandtl  number 
heat  transfer 

count  of  successive  solutions  of  the  momentum  equation 

in  the  procedure  of  solution 

radial  distance  from  axis  of  revolution 

radius  of  body  of  revolution 

recovery  factor,  defined  by  eq.(7.10) 

x  dr 
r  dx 


Reynolds  number  = 


o  u  x 
^e  e 


Stanton  number,  defined  by  eq. (6.120) 

transverse  curvature  term  defined  by  Eq. (6.32 )  except  in 

Sections  7.3  and  J.k  where  it  is  absolute  temperature 

x  component  of  velocity,  =  f ' 

e 

y  component  of  velocity 
vector  velocity,  eq.(A-5) 

distance  along  body  surface  measured  from  stagnation  point 
axial  or  chordwise  distance 
body  force  in  x-direction 
distance  normal  to  x 

nondimens ional  normal  distance,  eq.(6.106) 
body  force  in  y-direction 
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angle  between  normal  to  the  surface  y  and  xhe  radius  r 

Hartree's  measure  of  free-streara  velocity  distrioution  in 

similar  flow,  t3  =  — 

ratio  of  specific  heats,  Cp/cy 

thickness  of  boundary  layer 

displacement  thickness 

dimensionless  displacement  thickness,  eq. (6.103) 
specified  accuracy  in  program 
transformed  y  coordinate,  eq. (6.4) 

either  a  very  large  q,  that  is,  outside  the  boundary  layer 
or  q  at  the  edge  of  the  boundary  layer,  i.e.,  =  q 
momentum  thickness 

dimensionless  momentum  thickness,  eq. (6.110) 
dynamic  viscosity 
kinematic  viscosity 

variable  of  integration,  or  in  Section  7*9,  a  transverse 
curvature  parameter 

function  in  energy  equation  defined  by  eq. (6.50) 

density 

shear  stress 

transformed  stream  function  =  f  —  q,  see  eq.(6.j6) 

in  Section  6.1,  a  stream  function  defined  by  eqs.(6.5  and  6.7); 

in  all  other  Sections,  the  enthalpy  function  g  —  1,  eq.(6.45) 

SUBSCRIPTS 

evaluated  at  adiabatic  wall 

evaluated  at  outer  edge  of  boundary  layer 
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n 

rec 

ref 

stag 

t 


v 

00 
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Primes 


evaluated  at  station  n 

evaluated  at  recovery  temperature,  Trec  *  Tad 
evaluated  ar  reference  enthalpy 
evaluated  at  stagnation  point 
evaluated  at  total  temperature 
evaluated  at  vail 

evaluated  at  a  reference  condition,  or  see  t)n 

usually  means  evaluated  at  initial  condition,  x  «  0. 
One  exception  is  rQ 

denote  differentiation  with  respect  to  t|.  *  f',  etc. 
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5.0  INTRODUCTION 


The  present  work  was  started  because  there  were  many  problems  of  laminar 
boundary-layer  flow  for  which  no  satisfactory  solutions  had  been  found. 

Among  them  were  questions  on  high-speed  heat  transfer  and  the  effects  of  heat 
transfer  and  wall  mass  transfer  on  drag.  Answers  to  these  questions  require 
an  exact  general  method  for  solving  the  boundary-layer  equations.  Available 
methods  for  solving  the  boundary-layer  equations  appeared  inadequate.  Integral 
methods  can  treat  general  flows,  but  they  give  only  approximate  solutions. 
Similar-flow  methods  are  accurate  but  are  restricted  to  special  pressure  dis¬ 
tributions.  When  the  work  was  started,  finite-difference  methods  appeared  to 
require  unreasonably  long  computing  times  for  accuracy.  Therefore,  studies 
were  made  to  develop  a  practical  method  for  solving  "exactly1 1  the  complete 
equations  of  compressible  boundary-layer  flow  in  two -dimens ions.  The  sense 
of  "exactly"  as  used  here  is  that  the  solution  approaches  the  exact  as  the 
step  length  in  the  calculation  procedure  approaches  zero.  The  objective  of 
the  studies  was  to  find  a  method  capable  of  obtaining  solutions  for  arbitrary 
values  of  (l)  pressure  distribution.  (2)  wall  mass-transfer  distribution, 

(3)  gas  properties,  and  (k)  wall  temperature  distribution.  The  only  re¬ 
striction  on  the  flow  was  that  it  be  two-dimensional  or  axially  symmetric. 

The  problem  was  approached  by  first  finding  a  method  that  gave  an  ac¬ 
curate  solution  for  the  general  incompressible  boundary-layer  equation. 

An  accurate  and  rapid  method  was  developed  and  is  presented  in  References  1, 

2,  and  3.  The  method  is  a  modification  of  the  Hartree-Womersley  technique. 

It  consists  of  replacing  the  partial  derivatives  with  respect  to  the  flow 
direction  by  finite  differences,  while  re  taining  the  derivatives  in  a  direction 
normal  to  the  boundary,  so  that  the  partial  differential  equation  becomes 
approximated  by  an  ordinary  differential  equation.  The  method  was  programed 
on  an  electronic  computer,  and  solutions  for  a  wide  variety  of  flows  were 
calculated  and  are  presented  in  References  1  and  2.  The  large  number  of  cal¬ 
culations  and  their  comparison  with  other  exact  and  approximate  solutions 
establish  that  the  method  is  rapid,  highly  accurate,  and  powerful.  It  is 
capable  of  solving  any  flow  problem  for  which  the  incompressible  boundary- 
layer  equations  themselves  are  valid.  The  mathematics  of  the  solution,  that 
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is,  the  behavior  of  the  process  of  solution,  the  nature  of  solution,  its 
difficulties,  etc.,  are  reported  in  Reference  3.  Also,  the  reasons  for 
choosing  the  present  method  rather  than  the  more  conventional  finite-difference 
methods  are  discussed. 

After  the  method  of  solution  had  proved  successful  for  incompressible 
flow  it  was  extended  and  modified  to  solve  the  fully  general  laminar  boundary- 
layer  equations  for  compressible  flow.  The  only  restrictions  are  that  the 
flow  be  two-dimensional  or  axially  symmetric  and  that  no  dissociation  occur  . 

The  purpose  of  this  report  is  to  describe  this  method  for  solving  the 
compressible  flow  equations  and  establish  its  validity.  Section  6.0  develops 
the  equations,  by  starting  with  the  Navier-Stokes  equations.  The  second-order 
effect  of  transverse  curvature  is  retained.  Steps  in  the  procedure  for  solving 
the  equations  are  described,  and  the  reasons  for  choosing  them  over  other 
possibilities  are  discussed.  Section  J.O  consists  of  the  results  of  calcula¬ 
tions  of  a  wide  variety  of  flows.  They  establish  the  accuracy  of  the  method. 
Wherever  possible  comparisons  of  the  calculations  are  made  with  other  methods 
and  with  experiment. 
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6.0  DESCRIPTION  OP  METHOD  OF  SOLUTION 


6.1  Boundary-Layer  Equations 

The  problem  to  be  considered  is  axis yrame trie,  steady, equilibrium  flow 
about  a  body  of  revolution.  The  equations  necessary  to  describe  such  a  flow 
are  those  of  continuity,  momentum,  and  energy.  Also  the  equation  of  state 
and  relations  describing  the  fluid  properties  such  as  viscosity  and  specific 
heat  are  required. 

The  problem  will  be  restricted  to  high  Reynolds  number, so  that  the  Navier- 
Stokes  and  energy  equations  are  simplified  to  essentially  the  form  of  the 
boundary-layer  equations  as  originally  developed  by  Prandtl  in  190^.  The 
second-order  effect  of  transverse  curvature  that  wan  neglected  by  Prandtl  will 
be  considered  here.  Van  Dyke  considered  all  second-order  effects  identified 
as:  transverse  curvature ,  longitudinal  curvature,  slip,  temperature  jump, 
entropy  gradient,  stagnation  enthalpy  gradient  and  displacement  in  a  recent 
article  for  the  special  case  of  a  blunt  body  (Reference  b).  As  pointed  out 
by  Van  Dyke,  logically,  all  second-order  effects  should  be  considered  con¬ 
currently,  but  consideration  of  all  these  effects  is  beyond  the  scope  of  the 
present  work.  The  effect  of  transverse  curvature  is  included  because  of  its 
importance  in  predicting  boundary-layer  growth  on  long  slender  bodies,  such 
as  certain  missiles. 

The  basic  notation  and  scheme  of  coordinates  are  shown  in  figure  1, where 
u^  is  a  reference  velocity  and  ug(x)  is  the  velocity  of  the  main  flow  just 
outside  the  boundary  layer.  The  tern  He  represents  the  total  enthalpy  out¬ 
side  the  boundary  layer  and  is  constant.  Local  enthalpy  outside  the  boundary 
layer  hg  is  given  from  the  relation 

12  12 
He  =  he  +lue  =  h„  +  \  < 

The  coordinates  are  a  curvilinear  system  in  which  x  is  distance  along  the 
surface  measured  from  the  stagnation  point.  The  dimension  y  is  measured 
normal  to  the  surface.  Within  the  boundary  layer  the  velocity  components  are 
u  and  v,  being,  respectively,  in  the  x  and  y  directions.  The  body  radius 
is  rQ. 
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Figure  1-  Boundary  layer  on  a  body  of  revolution.  Coordinate  system. 


The  basic  steady- flow  equations  for  a  laminar  boundary  layer  in  the  above 
coordinate  system  are  developed  in  Appendix  A.  They  are,  including  the  trans¬ 
verse-curvature  terms: 

CONTINUITY 


r  ^  (r  p  U)  +  ^  (r  P  V) 


MOMENTUM 


KNERGY 


„  du  .  „  Su  1  dp  ^  u  Sr  Su  ^  5  /  Su  % 

r  sh  sh  i  i  Sr  r  u  sh  i,  Su  "I 

p  Lu  ST +  v  57  J  =  F  57  [  +  ^ l“p?}  u  57  J 

+  57|£  ^  ^(1“^  }  u37  j 


(6.1) 


(6.2) 


(6.3) 
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Equations  (6.2)  and  (6.3)  differ  from  the  forms  usually  obtained  when  the 
Prandtl  boundary-layer  approximations  are  made, because  they  contain  the  trans¬ 
verse  curvature  terms: 


and 


ji  dr  du 

r  5y  5y 


in  (6.2) 


JL  dr 
r  3y 


jfrSr *■*-£ »« 


du 

37 


in  (6.3) 


6.2  Transformed  Boundary-Layer  Equations 

Whereas  the  above  equations  could  perhaps  be  solved  by  the  method  to  be 
presented,  they  will  be  transformed  to  a  more  convenient  coordinate  system. 
Flilgge-Lotz  and  Blottner  have  solved  the  above  equations  essentially  as  they 
are  written,  though  simplified  relations  for  the  fluid  properties— viscosity 
and  Prandtl  number—  were  used  (Reference  6).  The  equations  as  written  have 
several  disadvantages;  for  example,  they  may  be  singular  at  x  =  0,  and  the 
boundary-layer  thickness  varies  greatly  with  distance. 

These  difficulties  can  be  overcome  by  stretching  of  the  coordinate 
normal  to  the  wall  y.  The  transformation  that  has  proven  most  effective 
for  accomplishing  the  stretching  is  that  of  Howarth-Dorodnitsyn,  where 

y 

%  =f  p  ay 
0 


It  is  felt  that  a  second  requirement  of  the  transformation  of  the  y-coordinate 
be  that  it  should  reduce  to  the  incompressible  transformation  that  was  studied 
extensively  by  the  authors  in  References  1,  2,  and  3«  Therefore,  consider  the 
transformation 


n  = 


(6.4) 
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Also,  it  Is  convenient  to  describe  a  stream  function  f  such  that 


a 


p 

P  r  v 


(6.5) 


Furthermore,  as  was  done  in  incompressible  flow,  it  is  convenient  to  introduce 
a  dimensionless  stream  function  f  such  that 


bf 


f* 


u 


(6.6) 


The  relation  between  f  and  f  is 


♦  -  Vw1*  f 


(6.7) 


In  order  to  transform  the  boundary-layer  equations  above  to  x,  T)  -  coordi¬ 
nates,  the  following  relations  are  used: 


*  (st)v  +  (&l  (M 


du 


a  i  n  a  .  i  n  e  a  , 


■  (fcl  (&),  * 

_  /  %  ’  3 

'  0  Sn 


A 


dy 


'  y 

Jpdy 


a 

sr 


(6.8) 


(6.9) 


The  continuity  equation  is  automatically  satisfied  by  the  definition  of 
the  stream  function. 

The  momentum  equation  (6.2)  becomes  the  following  in  the  transformed 
plane: 
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(6.io) 


where : 


primes  denote  differentiation  with  respect  to  tj, 
P  P 


ff" 


p  h 
He  e 


(6.11) 


c  = 


Poo  ^00 


“Oh 
Me  ^e 


(6.12) 


du 
x  e 


P  =  =  Pressure  gradient  parameter  (6.13) 


R  S  r"  §x  =  ^a^ius  parameter 


For  equation  (6.10)  the  boundary  conditions  are: 


(6.14) 


t)  =  o  : 


f«(0)  =  f;  =  0 

f  (0)  = 


(6.15) 


If  the  wall  is  impermeable,  fy=  0;  but  with  flow  through  the  surface 


d(*r) 

fc”t-TT 


Then 


<*r)w  = 


A 

-f  0v(*>  vw(6)  rw(5)  d* 


In  the  incompressible  reports  (References  1,  2,  and  3)  M  was  used  for  the 
pressure  gradient  parameter  but  it  is  changed  here  to  prevent  any  confusion 
with  Mach  number. 
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Note  that  v  positive  corresponds  to  blowing;  negative  to  suction.  The 
w 

outer  boundary  conditions  are 

ij  — ►  oo:  f 1  — *  1 

f“—*0 


A  few  of  the  more  important  properties  of  eq. (6.10)  are  noted.  If  the 

edge  velocity  is  of  the  form  u  =  c.xP,  (x/u  )(du  /dx)  is  identically  P.  If 
R  ei.ee 

r  =  c^x  ,  R  =  (x/r)(dr/dx).  If  P  and  R  are  constants  it  can  be  shown  that 
the  equation  is  independent  of  x  and  provides  the  so-called  similar  solutions. 


Equation  (6.10)  has  several  advantages  over  other  possible  forms.  Other 
forms  are  singular  at  x  =  0  and  require  that  an  initial  profile  be  specified, 
usually  some  distance  from  the  leading  edge;  but  in  (6.10)  the  term  containing 
the  x-derivatives  disappears  at  the  start  of  the  flow,  and  the  solution  can 
be  started  with  a  similar  flow.  Finally,  the  equation  reduces  to  the  incom¬ 
pressible  form  which  was  studied  extensively  in  References  1,  2,  and  3.  Thus, 
not  only  doe3  form  (6.10)  have  the  advantages  that  sure  discussed  in  these 
references  for  incompressible  flow,  but  also  the  large  number  of  flows  studied 
there  may  be  used  as  a  check  on  the  present  method. 


To  transform  the  energy  equation  (6.3)  to  the 
the  enthalpy  ratio  g 


and  let 


x,  ii  -plane,  first  define 
(6.17) 


(6.18) 


Then  (6.3)  becomes 
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1 

~r 


-fee' 


H 


C(1 


P  +  1 
2 


f  g’  +  x 


r 

I  f' 


(6.19) 


In  the  solution  of  (6.19)  either  the  wall  temperature  or  the  heat  trans¬ 
fer  at  the  wall,  corresponding,  respectively,  to  and  g^  will  be  specified. 
The  outer  boundary  condition  is 

g(il)—»l  as  r)  — ♦  oo  (6.20) 


6.3  Fluid  Properties 

Fluid  properties  that  appear  in  the  momentum  and  energy  equations  are 
density  p,  viscosity  u,  and  Prandtl  number  Pr.  These  equations  were  de¬ 
veloped  to  be  valid  so  long  as  the  fluid  is  in  equilibrium,  that  is,  the  fluid 
properties  are  functions  only  of  local  conditions  —  pressure  and  enthalpy. 
Previous  investigators  have  usually  used  simplified  laws  for  these  properties, 
such  as,  a  power  law  or  Sutherland's  law  for  viscosity  and  a  constant  value 
of  the  Prandtl  number.  Whereas  it  is  known  that  these  laws  do  not  accurately 
represent  the  properties  over  the  entire  flight  regime  of  interest,  the  effect 
of  these  inaccuracies  on  the  solution  of  the  boundary-layer  equations  has  never 
been  investigated.  The  present  method  of  solution  has  been  developed  so  that 
arbitrary  fluid  properties  may  be  used.  That  is,  they  are  inputs  in  the 
program  in  the  form  of  formulas  or  tables  as  functions  of  local  enthalpy  and 
pressure.  There  is  no  restriction  on  the  fluid  to  be  considered,  that  is, 
the  flow  may  be  in  either  air,  helium,  water,  or  some  other  medium.  All  that 
is  required  is  that  the  fluid  be  in  equilibrium  and  its  properties  be  specified 
in  the  proper  form. 

All  of  the  flows  to  be  presented  in  this  report  are  for  air  without  dis¬ 
sociation.  Relations  for  fluid  properties  that  will  be  vised  have  been  chosen 
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from  Reference  7.  Reasons  for  this  choice  are:  (l)  The  relations  correlate 
well  with  the  known  properties  of  air  throughout  the  sensible  flight  regime, 
up  to  velocities  of  29,000  ft/sec  and  an  altitude  of  250,000  ft.  (2)  The 
relations  are  presented  as  functions  of  local  enthalpy  and  pressure,  whereas 
most  other  reports  present  the  properties  as  function  of  temperature  and 
pressure.  (5)  The  reference  also  considers  equilibrium  dissociating  air  and 
the  relations  presented  are  compatible  when  either  dissociation  or  no  dis¬ 
sociation  occurs.  Whereas  the  present  report  is  restricted  to  nondissociating 
gases,  it  is  planned  to  extend  the  method  presented  here  to  dissociating  gases 
in  equilibrium.  Thus  the  relations  of  reference  7  would  be  consistent  for 
such  an  extension  to  flow  in  air. 


Cohen  (Reference  7)  presents  formulas  for  density  and  the  product  of 
density  and  viscosity  and  a  table  of  Prandtl  numbers  all  as  functions  of  a 
local  enthalpy  ratio  h/href  and  static  pressure  p.  He  chose  href  to  be 
about  the  average  value  of  total  enthalpy  encountered  in  the  sensible  flight 
regime 

href  *  2.119  x  108  ft2/sec2  (6.21) 


The  enthalpy  ratio  can  be  obtained  from  the  functions  of  f  and  g  of  the 
momentum  and  energy  equations  from  the  relations: 


H 


2  h 


'ref 


(6.22) 


He 

=  STT  _ 

ref 


-  f  »2 


u 


u 


.2  2h. 


ref 


(6.25) 


2  2  2 
Mote  that  ^e/\ef  Bad  u{o/href  are  constants  for  the  flow,  ug/um  is  a 

function  of  x  only,  and  g  and  f  *  are  functions  of  both  q  and  x. 


Cohen's  formulas  for  the  variation  of  density  with  enthalpy  in  the 
boundary  layer  (recall  that  there  is  no  pressure  change  in  the  boundary  layer) 
can  be  written  as 
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(6.24) 


Pe  (h/hyef )°*6123  “  0.0455283 
p  (he/href  )0'6123~  °* ^55283 

With  regard  to  the  accuracy  of  (6. 24), Cohen  states  1 ’the  function  fits  the 
data  reasonably  well  over  the  range  0.0152  i  h/href  S  2.0“  for  pressures 
in  the  range  10“  to  10  atmospheres.  “The  maximum  deviation  in  this  range 
is  about  t  25  percent  at  low  enthalpy  and  the  average  deviation  for  all  data 
is  about  ±  5  percent.  Agreement  is  best  in  the  range  0.2  <  h/href  <  1.6“ . 

Viscosity  appears  in  the  boundary-layer  equations  as  the  product  of 
viscosity  and  density.  This  product  is  a  function  of  enthalpy  only  across 
the  boundary  layer  ,  and  its  ratio  C  according  to  Reference  7  is 


p  it  (hg/h ref)0*5329  “  0.020856 

pe^e  (h/href)0*^^  -  0. 020856 


(6.25) 


The  variation  of  C  outside  the  boundary  layer  is  also  required  and  is  given 
by 


where 


P.A, 

Pe^e 


(Vhref)0’5529  ”  °*020856  I  Pm  \°*"2 

(h 09/href  -  0.020856  \  pe  / 

(6.26) 


H 


u 


‘ref 


2h 


ref 


With  regard  to  the  accuracy  with  which  (6.25)  and  (6.26)  agree  with  data, 
Cohen  states ; ' 1  Maximum  deviation  is  about  ±  8  percent  and  the  average  de¬ 
viation  in  the  entire  range  0.0152  s  h/h^ef  S  2.0  is  about  i  3  percent.” 
For  flows  where  W\ef  is  1®SB  than  0.015,  a  linear  viscosity  law  is  used 
and  the  density  is  taken  as  inversely  proportional  to  local  enthalpy. 


Reference  7  presents  the  variation  of  Prandtl  number  with  enthalpy  in 
table  form  (Table  1,  page  26  of  Ref. 7).  Whereas  it  is  possible  to  use  this 
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table  as  an  input  in  the  computer  program  for  solving  the  boundary-layer 
equations ,  a  simpler  and  more  accurate  method  is  to  fit  relations  to  the 
tabular  values.  It  has  been  found  that  the  following  three  least-square 
relations  fit  the  tabular  values  with  an  error  of  less  than  0.22  percent 
over  the  range  of  data. 


For  0.005 


where 


For  0.075 


where 


i  h/href  S  0.075, 

6 

Pr  *  0.77  +  V  1 

i  (-r-  —  0.005 )n 

n  ref 

a1  -  —  6.18253 

a^  =  -  642822.0 

ag  -  -  147.92^5 

a5  =  +  8.00559  x  10( 

aj  =  +  21609.81 

ag  =  —  3.66200  x  10 

-  h/href  S  °-30' 

6 

pr  =  0.7374  +y 

Vh/hr«f  ~  0-°75)n 

n=] 

L 

bx  =  +  2.009 

\  -  “  5319.69 

b2  =  —  45.112 

bj.  =  +  10613.04 

b?  =  +  524.907 

bg  =  —  13410.82 

(6.27a) 


(6.27b) 


For  0.50  S  5  2.0, 


Pr  =  0 


2 

•  755  cn( 


(h/href  “  °-50)n 


(6.27c) 


where 


c]_  =  -  0.1299 
c2  =  +  0.05757 
Cj  =  +  0.001323 
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Figure  2- Variation  of  Prandtl  number  with  enthalpy. 


These  relations  and  the  tabular  values  from  Reference  7  are  plotted  in  figure  2. 
At  very  low  enthalpies,  it  is  seen  that  the  Prandtl  number,  both  Cohen's  data 
and  the  fitted  curves,  increases  rapidly  as  enthalpy  decreases.  It  is  believed 
that  Prandtl  number  should  have  a  value  of  0.72  at  the  lower  enthalpies. 
Therefore,  if  <  0.015,  this  value  is  vised  in  the  program. 

Boundary  layer  in  flows  of  fluids  other  than  air  or  in  air  with  other 
fluid-property  laws  can  be  handled  by  replacing  the  relations  of  this  section 
by  ones  appropriate  to  the  flow  being  studied. 

6.4  Choice  of  Procedure  for  Solution  of  Boundary-Layer  Equations 

The  method  of  sol"tion  will  be  similar  to  that  used  in  the  study  of  the 
incompressible  bourn*  _ayer  (References  2  and  3).  The  x-derivatives  in  the 
momentum  and  energy  equations  are  replaced  by  finite  differences, so  that  the 
partial  differential  equations  are  approximated  by  ordinary  differential 
equations .  Then  the  problem  of  solution  is  essentially  to  find  the  unknown 
boundary  conditions  at  the  wall  that  satisfy  the  known  outer  boundary  con¬ 
ditions.  This  is  done  by  a  cut-and-try  procedure,  which  is  described  in 
Sections  6.5  and  6.6.  The  momentum  equation  (6.10 )  and  the  energy  equation 
(6.19)  are  interdependent  and  must  be  solved  simultaneously. 

Several  procedures  for  solving  the  equations  simultaneously  are  possible. 
Three  are: 

I.  Starting  with  assumed  boundary  conditions  at  the  wall,  solve  the 
two  equations  simultaneously  with  the  appropriate  fluid  properties. 

A  cut-and-try  procedure  would  be  used  on  the  wall  values  until  the 
outer  boundary  conditions  are  satisfied. 

II.  To  get  started,  assume  an  enthalpy  distribution  (and  thus  the  fluid 
properties )  and  calculate  the  momentum  equation.  Again  a  cut-and- 
try  procedure  would  be  required.  Values  of  the  stream  function  f 
and  the  velocity  f '  from  the  first  solution  of  the  momentum  equa¬ 
tion  would  be  used  in  the  solution  of  the  energy  equation  and  .the 
fluid  properties.  The  new  fluid  properties  would  be  used  again  to 
solve  the  momentum  equation.  This  iterative  procedure  would  be  con¬ 
tinued  until  convergence  of  the  solution  is  obtained.  A  diagram  of 
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the  procedure  is  sketched  below 


III.  Again, to  get  started,  the  fluid  properties  are  assumed  to  get  a 
first  solution  of  the  momentum  equation.  The  energy  equation  is 
then  calculated,  but  by  using  the  assumed  fluid  properties.  By 
use  of  this  assumption,  the  energy  equation  is  linear,  and  thus 
its  solution  is  considerably  simplified.  Then  the  calculated 
enthalpies  are  vised  to  determine  new  fluid  properties, which  are 
used  to  solve  anew  the  energy  equation.  Once  convergence  of  the 
fluid  properties  is  obtained,  the  momentum  equation  is  solved  a 
second  time.  This  double-iteration  procedure  is  continued  until 
convergence  of  both  the  momentum  and  energy  equations  is  obtained. 
The  procedure  is  diagramed  on  the  next  page. 

Studies  indicate  that  the  solution  will  be  simpler  and  probably  more  accurate 
if  derivatives  of  the  fluid  properties,  such  as  dC/dtj,  d(pe/p)/dn,  etc.,  do 
not  have  to  be  evaluated  in  the  solution.  This  forces  the  use  of  Procedure 
III  above,  for  consider  the  energy  equation  (6.19) 


26 


First 

Solution 

of 

Momentum 

Equation 


Second 

Solution 

of 

Momentum 

Equation 
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r~  HI.  PROCEDURE  FOR  SOLVING  EQUATIONS  SIMULTANEOUSLY 


< 


Iterative  procedure  is  con¬ 
tinued  until  calculated  va¬ 
lues  of  enthalpies  converge. 


Above  iterative  procedure 
repeated  until  calculated 
values  in  momentum  equation 
converge. 
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J 

L  2  J 
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(6.19) 


from  using  values  of  dc/dt)  one  must  solve 
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Pr 
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s'  +  ~TT  )  ft 
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fP  +  1 


+  R 


(6.28) 


t  g 


dt)  +  K 
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Thus  the  left-hand  side  of  (6.28)  can  be  determined,  but  the  values  of  C  and 
Pr  depend  on  g,  which  must  be  determined  by  Integrating  g*.  An  Iterative 
procedure  could  be  used  to  determine  the  value  of  C,  Pr  and  g‘,  but  this  can 
be  avoided  by  using  Procedure  III. 

Detailed  steps  in  the  procedure  of  solution  are  given  in  Section  6.9 
and  Appendix  B . 


6.5  Method  of  Solution  of  Momentum  Equation 


Before  attempting  its  solution,  equation  (6.10)  will  be  rewritten  for 
simplification.  First, consider  the  first  term  in  (6. 10) 


-r- 1  !r(0rf”) 


The  radius  r  is  a  function  of  both  x  and  t),  being  defined  by 


r  =  r  +  y  cos  a  =  r  +  cos  a 
o  o 


<6-29> 


where  rQ  is  the  local  radius  of  the  axisymmetric  body 

a  is  the  angle  between  the  normal  to  the  surface  and  the  radius  r 
(see  Fig.  Al). 

It  would  be  difficult  to  handle  the  integral  in  (6.29)  witliin  the  differential 
in  equation  (6.10).  Therefore  rewrite  the  first  term  of  (6.10)  as 


1  1  a 

IT  r  Sn 
00 


(C  r  f") 


l_  8 

c 

00 


(c  f")  +  -§-  f"  7 

00 


Sr 


(6.50) 


The  second  tern  on  the  right-hand  side  is  the  transverse- curvature  term  and  is 
negligible  when  the  boundary -layer  thickness  is  small  compared  to  the  body  radi¬ 
us.  Thus  there  is  a  second  advantage  to  writing  the  first  term  of  (6.10)  in  two 
parts, in  that  the  transverse- curvature  tern  stands  by  itself.  For  further 
simplification  introduce  T,  a  transverse  curvature  parameter,  defined  such 


that 
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Substitution  of  the  value  of  r  given  in  (6.29)  gives 


Now  T  may  be  written  as 


(6.31) 


(6.32) 


The  radius  parameter  R  may  be  written  as 

R  .  x  dr  [  x _ [  ^o  d  cos  a  "I 

*  F  dx  ”  rQ  +  y  cos  a  dx  y  dx  J 


(6.33) 


Now  if  the  boundary  layer  is  small  with  respect  to  the  body  radius,  as  is 
usually  assumed  in  boundary-layer  flow,  R  is  simply 
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dr 

o 

£T 


(6.33a) 


R 


X 


But  if  it  is  not,  the  y  terns  in  (6.33)  may  affect  the  solution  of  the 
momentum  equation.  The  last  tern  in  (6.33)  is 


dr  , 

o  ,  „  d  cos  a 

dx  y  dx 


dr 

^T-y  sin  a 


da 

33c 


da 

"y  Hx  _ 


since  sin  a  =  dro/dx.  But  da/dx  is  the  reciprocal  of  the  local  longitudinal 
radius  of  curvature,  and  this  radius  must  he  large  with  respect  to  y  if 
the  boundary-layer  approximation  used  in  developing  the  momentum  and  energy 
equations  is  to  hold.  That  is,  if  the  local  radius  of  curvature  is  small 
with  respect  to  y,  the  second-order  effect  of  longitudinal  curvature  becomes 
of  importance,  and  the  approximation  of  no  pressure  change  across  the  boundary 
layer  is  no  longer  true.  Therefore  y(da/dx)  <  <  1  and  (6.33)  may  be  written 
as 
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x_  dr 
r  dx 


H3T 


1  +  ^  cos  a 
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(6.33b) 


Also  for  simplicity  the  symbol  N  is  introduced,  defined  as 


N 


P  +  1 
2 


+  R 


(6.34) 


The  momentum  equation  can  now  be  written 
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Ir  (° f")  -  - T  r 0  f"  -  c.r<T- f,2) "  °." rf" +  c-x  [ft  s’  - f”  i] 

(6.35) 

It  will  be  solved  In  the  same  manner  as  the  Incompressible  equation  was  in 
References  2  and  3«  Values  of  the  fluid  properties  C  and  p  /p ,  which  of 
course  were  not  needed  in  the  incompressible  problem,  are  assumed  here  to  be 
given  by  the  most  recent  solution  of  the  energy  equation.  For  the  very  first 
solution  they  are  assumed. 


In  the  authors'  studies  of  the  incompressible  boundary  layer  it  was  found 
that  round-off  errors  in  the  computer  program  could  be  reduced  by  substitution 
of  cp 1  =  f  *  —  1  into  the  moment™  equation  (Reference  3).  The  same  substitu¬ 
tion  will  be  made  here.  The  reason  that  round-off  errors  are  reduced  by  the 
substitution  is  that  in  equation  (6. 35)  all  terns  approach  zero  as  t)  ap- 
proaches  oo ;  both  pg/p  and  f  approach  unity  and  the  round-off  error 
is  primarily  introduced  when  taking  their  difference.  The  substitution  is 

<p  =  f  -  n 

cp  *  *  f '  —  1 

<p»»  s  f»» 

cpitt  *  f».. 

Introduction  into  (6.35)  gives 
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(6.38b) 


cp*  * — ►  o 


The  region  of  solution  will  be  divided  into  x-wise  stations  as  shown  in 
figure  3«  At  each  station  the  x-derivatives  will  be  replaced  by  finite 
differences.  The  finite-difference  representation  is  described  in  Section 
6.7.  This  replacement  of  the  partial  derivative  with  respect  to  x  by  finite 
differences  results  in  an  ordinary  differential  equation  at  each  x-station  for 
the  momentum  equation.  Each  equation  must  be  solved  step  by  step  as  the  cal¬ 
culation  proceeds  in  the  x-direction.  The  equation  is  third  order  and  non¬ 
linear.  Solution  of  the  equation  is  made  difficult  by  both  the  nonlinearity 
and  the  fact  that  one  boundary  is  at  tj  =  00.  Ample  work  in  the  past  has 
proved  the  existence  of  a  solution)  therefore  it  is  sufficient  to  search  for 
the  correct  solution.  A  positive  method  for  doing  this  is  to  solve  (6.37)  as 
an  initial  value  problem  using  arbitrary  values  of  cp^'  as  a  third  boundary 
condition.  It  is  then  necessary  to  search  through  the  possible  values  of  cp^' 
to  find  the  one  that  satisfies  the  outer  boundary  condition  —  that  cp*  ap¬ 
proaches  zero  asymptotically  as  approaches  infinity.  The  procedure  for 
performing  the  search  is  described  below  but  first  consider  the  solution  of 
tne  equation  as  an  initial  value  problem.  First  determine 


n 

C  9M  =/^f  (C  9,,)  dT1  +  Cv  9w'  (6*59) 

0 

where  ^  (C  cp")  is  given  by  the  right  hand  side  of  (6.37).  is  known  or 

assumed  but  cp^ '  is  not.  It  will  be  found  by  the  searching  procedure  described 
in  the  following  paragraph.  The  method  of  integration  to  be  vised  in  (6.39)  is 
described  in  Section  6.8.  Other  quantities  needed  in  the  solution  of  (6. 39) 
are  given  by 


cp*  * 


(c  cp”) 
C 


cp* 


dii  -  1 


(6.40) 

(6.41) 
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(6.42) 


n 

cp  -  JV  dn  +  <pw 

o  * 

Briefly,  the  procedure  for  searching  for  the  correct  value  of  is  to 
first  try  values  of  until  the  solutions  for  cp*  are  bounded  in  a  specified 
region.  Then  an  interpolation  method  is  used  to  obtain  the  value  of  cp^'  that 
satisfies  the  outer  boundary  condition.  The  steps  in  the  procedure  are: 

1.  Try  <p»*  *  (q" ) .  where  the  latter  is  an  input  to  the  program 

v  w  input 

or,  more  conveniently,  is  the  value  at  the  previous  station.  Compute 

outward  to  determine  if  the  trial  solution  exceeds  cp'  *  0  or  not. 

Because  of  the  transformation  used,  the  value  of  q^*  will  generally 

remain  between  0  and  4.  A  value  of  cp"  =  0  corresponds  to 

w 

separation. 

2.  If  cp'  exceeds  0,  the  trial  value  of  qj’  is  high  and  a  second 
solution  is  computed  with  a  reduced  q^'.  This  procedure  is  continued 
until  both  a  high  and  a  low  value  of  cp£'  are  known. 

3-  Once  both  a  high  and  a  low  value  of  q^ '  are  known,  the  bounds  on  the 
correct  value  of  q^*  can  be  further  narrowed  by  splitting  the 
difference  between  the  upper  and  lower  bounds,  and  computing  again. 

4.  This  splitting  the  difference  can  be  continued  until  q^'  for  a  high 
solution  and  q^'  for  a  low  solution  agree  to  a  specified  number  of 
decimal  places.  This  procedure  is  positive  but  costly  in  computing 
time.  Studies  have  shown  that  the  searching  procedure  can  be  speeded 
up  considerably  with  no  loss  in  accuracy  by  the  following:  the  split- 
ting-the-difference  procedure  is  used  until  three  solutions  are  ob¬ 
tained  such  that  q1  at  is  between  the  bounds  of 

—  K  S  q'  (t^)  i  K.  (Both  and  K  are  inputs  in  the  computer 
program. )  At  least  one  of  the  three  solutions  must  be  high  and  one 
low.  A  three-point-interpolation  procedure  is  then  used  to  determine 
the  solution  that  satisfies  the  outer  boundary  condition  ?'(%,)  =  0. 
This  is  the  same  procedure  that  was  used  in  the  solution  of  the  in- 
compressible-boundary-layer  equations  and  reported  on  in  Reference  3» 

The  interpolation  procedure  is  as  follows.  Consider  the  typical  set  of 
trial  runs  shown  in  figure  4.  The  runs  are  for  incompressible  stagnation-point 
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Figure  4.-  Trial  solutions  for  various  values  of  <p«  •  Stagnation-point  flow,  P  =  1.0,  R  =  0. 


flow.  The  tries  were: 


Try 

<P" 

Nature 

Tries  used  for 

T)max 

W 

interpolation 

1st 

1.0 

Low 

2.50 

2nd 

1.5 

High 

1.95 

3rd 

1.25 

High 

1st 

5.0 

4th 

1.125 

Low 

5.55 

5th 

1.1875 

Low 

2nd 

5.0 

6th 

1.21875 

Low 

3rd 

5.0 

Interpolated 

1.232587 

The  program  first  tried  <p^'=  1*0,  which  it  found  to  he  low.  The  second  try 
of  1.5  was  high.  It  then  proceeded  to  split  the  difference  between  the  last 
high  and  low  tries  until  it  had  three  solutions  that  extended  all  the  way  to 
within  the  bounds  of  cp* )  ±  1.  These  three  solutions  were  the  3rd, 
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5th,  and  6th  tries.  Denote  them  as 


1st  solution 

">1 

<pM 

(c  cp"){ 

2nd  solution 

*2 

*2 

<P£* 

(C  q>M)£ 

3rd  solution 

■e 

VJJ 

n 

(pit 

(C  cp")$ 

Lagrangian  three-point  interpolation  is  used  to  determine  the  solution  which 

meets  the  outer  boundary  condition  ©’(▼)  )  =  0.  The  interpolated  solution  is 

00 

given  by  \ 


q>  (ti)  -  ^(n)  +  Ag  <p2(n)  +  a5  cp5(n) 
9»(tj)  =  ^  q>£<n)  +  Ag  cp^(rj)  +  A?  cp»(n) 


(6.43) 


and  a  similar  relation  for  <p*  •  (tj )  and  (C  cp")'  where  the  coefficients  are 
given  by 


<p£(n„)  ^(nco) 

*2  =  T«£(tJ  -  -  ^<03  f 

3 c  ' 


(6.44) 


The  solution  can  be  made  as  accurate  as  desired  by  restricting  the  values 
of  the  bounds  K.  Effect  of  K  on  the  accuracy  of  solution  is  discussed  in 
Section  7«2.  Typically  K  would  have  a  value  of  1  for  five-place  accuracy. 
Computing  time  required  to  obtain  this  accuracy  with  the  interpolation  method 
is  about  one  fourth  the  time  required  if  just  splitting  the  difference  between 
the  high  and  low  tries  were  used.  A  study  showed  the  three- point  form  of 
interpolation  to  be  considerably  more  accurate  than  the  two-point  form,  but 
no  great  gains  in  accuracy  were  obtained  by  using  more  than  three  points. 
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6.6  Method  of  Solution  of  Ihergy  filiation 


Before  solving  the  energy  equation  (6.19)  It  Is  convenient  to  substitute 
the  functions  T  and  cp  Introduced  In  the  preceding  Section  and  to  also 
Introduce  the  function  \|r  defined  as 


♦  -  g  -1 


(6.4 5) 


for  the  same  reasons  that  cp  was  Introduced  In  the  momentum  equation.  Substi¬ 
tution  of  the  functions  in  (6.19)  gives 

_  r  u2  . 

It  -w*1  +  ^r  cd-& )(*'  + i)<p"  s 


=  -t(-£  ) 

p 


-pj-  ♦'  +  -g£-  C(l  )(q>*  +  1)  cp" 


—  C  N(cp  +  Tj)t'  +  c  X 


(6.46) 


The  wall  boundary  condition  are: 
at  q  =  Oj 

*V  -  e»  - 1 

A  more  general  boundary  condition  could  be  defined  by 


(6.47) 


n  =  gw  +  kgyj 

but  usually  either  or  gj,  corresponding  respectively  to  wall  tempera¬ 
ture  and  wall  heat  transfer,  would  be  known. 


The  third  boundary  condition  is,  as  T)  — ►  oo; 

i)r  — ♦  0 
t'— ►  0 


(6.48) 
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The  method  of  solution  of  (6.46)  Is  similar  to  that  of  the  momentum 
equation.  The  region  of  solution  is  divided  Into  x-wise  stations  as  indicated 
in  figure  J.  Again  the  x-derivatives  are  replaced  by  finite  differences , which 
are  defined  in  the  next  section  of  this  report.  When  solving  (6.46)  values 
of  cp  and  its  derivatives  that  were  determined  from  the  previous  solution 
of  the  momentum  equation  and  fluid  properties  that  were  determined  from  the 
previous  solution  of  the  energy  equation  are  vised.  Once  a  solution  has  been 
obtained  the  fluid  properties  are  re-evaluated  and  the  energy  equation  is 
solved  again.  This  iterative  procedure  is  continued  until  convergence  of 
the  solution  is  obtained.  Details  of  the  iterative  procedures  for  solving 
both  the  momentum  and  energy  equations  are  given  in  Section  6.9.  With  these 
procedures  i]r  and  \|r*  are  the  only  unknowns  in  the  solution  of  (6.46)  and 
the  equation  is  linear. 


The  solution  of  (6.46)  is  as  follows.  Rewrite  it  as 
P, 


(tt)  s  tt'  =  _  t  ?r  -  CM^cp  +  tj  )i|f  *  +  C 
where  for  simplicity 


(6. 


0 


Pr 


r  +  —■  c(l  “ ]5jr  )(<?’  +  1)  9 


(6.49) 


(6.50) 


Equation  (6.49)  is  integrated  to  determine  rr  : 

n 


t  =  J  ir' 

0 


drj  +  TT 


(6.51) 


The  method  of  integration  is  the  same  as  for  the  momentum  equation  and 
is  described  in  Section  6.6.  From  (6.50) 

2 

r-  U  H 

v  =  — ir  c(i--Kr)(9'  +  1)  c?” !  (6.52) 


Pr 


which  may  be  integrated  to  determine  \j»  : 

4 

V  =  /V  dn  +  * 


(6.53) 
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Since  equation  (6.1+9)  is  treated  as  linear,  its  solutions  may  be  linearly 
combined.  It  will  be  solved  twice, and  the  two  solutions  combined  to  meet  the 
outer  boundary  condition.  The  exact  procedure  will  depend  on  whether  or 
gj  is  known. 

Case  1:  g^  is  known.  Both  solutions  begin  with  the  same  value  of 
the  one  Imposed  by  the  boundary  condition 


♦(0)  -  *w  =  g„  -  1 

First  equation  (6.1+9)  is  solved  by  vising  a  trial  value  of  ^  . 

The  solution  is  denoted  as 

♦jOi) 

If  is  sreater  than  zero, a  lower  value  of  f '  is  tried;  if  it 

is  less,  a  higher  value  of  .  The  second  try  is  denoted  as 

♦2^) 

The  two  solutions  can  be  added  to  produce  the  most  general  solution, 
which  can  be  made  to  meet  the  boundary  conditions.  The  general  solution 
is 

♦(1)  =  a  ^(n)  +  b  t2(n)  (6.51+) 

The  outer  boundary  conditions  are 


*(nj  =  A  ^(nj  +  b  t2(iw)  *  0 
♦'(To,)  =.a  ♦{(»!»)  +  B  ^(n,,)  =  0 


(6.55) 


and  also 


But 


so  that 


At  +  B 


A  +  B  =  1 
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Equation  (6. 55)  gives 


-»2(tO 

and  the  correct  solution  is  for  all  tj*s 

♦(n)  =  a  ^(n)  +  (1  -  a)  i|t2(t)) 

and  also 

t’(n)=  a  ¥{(n)  +  (1  -  a) 


(6.56) 


(6.57a) 

(6.57b) 


Case  2:  g^J  is  known. 

The  procedure  is  similar  to  Case  1,  but  now  the  energy  equation  is 
solved  with  two  trial  values  of  instead  of  g*  .  Again,  the  two 

trial  solutions  are  denoted  as  i|u(t))  and  ^(tj).  Relations  (6.56), 
(6.57a), and  (6.57b)  then  may  be  used  to  give  the  correct  solution. 


6.7  Finite-Difference  Representation  of  x-Derivatlves . 


The  fundamental  idea  for  the  method  of  solution  —  that  of  replacing  the 
x-derivatives  by  finite  differences  to  approximate  the  partial  differential 
equation  by  an  ordinary  differential  equation  —  was  advanced  by  Hartree  and 
Wcmersley  (Reference  8).  The  idea  was  applied  to  the  incompressible  boundary 
layer  by  Hartree  (References  9  and  10 )  and  the  authors  in  References  1,  2, 
and  3. 


Two  treatments  within  the  scheme  are  possible.  One  is  to  deal  in  terms 
of  the  differential  equation  at  a  point;  in  particular,  the  x-derivatives  at 
a  point  are  replaned  by  their  finite>-difference  equivalent*.  Another  treat¬ 
ment  is  to  deal  in  terms  of  mean  values  of  the  variable  9  or  ^  for  a 
region  of  finite  extent.  Both  methods  are  described  in  Reference  3.  Within 
both  of  the  methods  there  are  many  possible  representations  of  a  derivative 
by  finite  differences,  for  example,  two,  three,  or  more  points  may  be  used. 
Note  that  all  of  the  x-derivatives  that  appear  in  the  momentum  and  energy 
equations  are  only  of  first  order. 
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In  reference  3  the  authors  made  an  extensive  study  of  the  use  of  two- 
point,  three-point,  and  four-point  finite  differences  in  solving  the  in¬ 
compressible  -boundary-  layer  equations.  Both  the  concept  of  the  point  form 
and  of  the  mean  form  were  studied.  The  investigation  showed  that  whereas 
the  two-point  mean  representation  was  more  accurate  than  the  two-point  point 
form,  the  use  of  three  and  four  points  proved  the  point  forms  to  be  more 
accurate.  Solutions  obtained  by  the  mean  method  with  the  higher  number  of 
points  diverged  wildly  as  the  step  length  became  small.  The  same  type  of 
divergence  appeared  in  Reference  1  where  the  two-point  mean  form  was  used, 
but  occurred  there  at  much  smaller  step  sizes  in  x.  The  investigation  also 
showed  that  whereas  the  use  of  three  points  in  the  point  form  of  finite 
differences  gave  a  much  more  accurate  solution  than  the  use  of  two  points, 
further  increase  of  the  number  of  points  to  four  gave  no  great  increase  in 
accuracy.  As  a  result  of  the  investigation,  the  three-point  point  form  of 
representing  the  derivatives  was  chosen  in  the  present  method.  In  addition 
to  being  freer  of  oscillations  as  the  step  size  becomes  very  small,  the  point 
form  has  the  advantages  of  being  simpler  conceptually  and  simpler  to  program 
on  the  computer. 

The  basic  scheme  of  the  finite-difference  representation  is  diagramed 
in  figure  5*  The  space  is  divided  into  a  number  of  regions  bounded  by  lines 
Xjj,  xn_i>  Xn-2*  xn-3"  sPac*nS  A  x  need  not  be  constant.  Because  the 

momentum  equation  and  energy  equation  are  parabolic  in  x,  the  problem  must 
be  solved  by  proceeding  in  the  direction  of  positive  x.  It  is  assumed  that 
the  solution  has  been  found  at  all  previous  stations  up  to  and  including 
xn_1,  which  of  course  means  that  cp(Tj)  and  ^(tj)  and  their  derivatives  are 
fully  known  at  these  stations.  The  quantity  cp' (tj)  typically  has  the  appear¬ 
ance  shown  in  figure  5*  The  problem  is  to  find  the  solution  q>(  “n )  and  \Jr ( rj ) 
at  the  new  station  xn. 

Whereas  usually  three-point  finite  differences  will  be  used,  at  the  start 
of  a  solution  only  two  points  are  available  and  the  two-point  form  must  be 
used.  Also,  when  there  are  discontinuities  in  the  boundary  conditions,  it  may 
be  more  accurate  to  use  the  two-point  form.  The  finite  difference  repre¬ 
sentations  are: 
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The  error  here  is 


A 


(*n  ~  xn-l^xn  ~  Xn-2^  d3  cp(j) 

5  ax5 

where  here  |  is  some  value  of  x  in  the  interval  xQ  —  ^n_2  • 

In  solving  the  boundary- layer  equations  the  quantity  dcp/dx  in  both  the 
momentum  equation  (6.37)  and  the  energy  equation  (6.46)  is  replaced  by  either 
(6.58)  or  (6.59).  The  other  x-derivatives  8cp'/^x  &ud  &i|r  /dx  are  replaced 
by  similar  expressions.  When  these  substitutions  are  made,  it  is  assumed 
that  all  of  the  other  quantities  in  the  boundary-layer  equations  are  evalu¬ 
ated  at  xq.  The  equations  are  then  ordinary  differential  equations  in  t) 
with  the  variable  quantities  cp,  c?1,  and  i|r  at  the  n-1  and  n-2  stations. 
Step  length  Ax  is  not  a  primary  parameter*  instead,  x/a  x  is. 

In  solving  the  boundary-layer  equations  the  calculation  must  start  at 
x  =  0.  For  the  x  =  0  station  the  terms  with  x-derivrAives  in  both  the 
momentum  and  energy  equations  disappear.  At  the  second  station  the  two-point 
form  of  the  finite  differences  is  used, but, at  all  stations  farther  downstream 
the  three-point  form  may  be  vised.  The  error  in  the  three-point  form  is  like 

iA_ 

5  dx5 

as  compared  to 

A  x 

T- 

for  the  two-point  form.  Therefore  in  order  to  have  the  same  accuracy  in  the 
solution  at  all  stations  the  step  size  at  the  second  station  must  be  suitably- 
reduced.  In  practice  these  errors  near  the  leading  edge  will  probably  be 
small  due  to  the  fact  that  the  flow  is  nearly  similar  there,  and  thus  the 
x-derivatives  are  small.  But  in  any  case  A  x  can  be  kept  small  here,  in 
order  to  keep  accuracy  high. 

For  further  discussion  of  application  of  Har tree - Womer s ley's  method  for 
solution  of  the  boundary-layer  equations  and  the  associated  errors,  see 
References  1  and  3. 
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6.8  Method  of  Integration 


The  overall  method  of  solution  of  both  the  momentum  and  energy  equations 
is  outlined  in  Sections  6.5  and  6.6.  The  problem  of  solution  is  essentially 
one  of  integration.  There  are  several  methods  for  performing  the  integration 
that  are  available  as  subroutines  on  the  computer,  for  example,  Milne's 
fourth-order  predictor-corrector  method,  but  because  of  their  generality, they 
require  long  computing  times  to  solve  the  present  problem.  Therefore  the 
authors  made  a  study  of  various  techniques  that  were  available  for  performing 
the  integration  as  applied  to  the  incompressible  problem.  The  results  of 
this  study  are  reported  in  Reference  3*  A  method  was  developed  there  that 
appears  to  be  both  highly  accurate  and  rapid.  It  uses  the  Falkner  extra¬ 
polation  formulas  and  the  Adams  interpolation  formulas  which  are  described 
on  pages  116-131  of  Reference  11.  Furthermore, as  a  result  of  the  study,  the 
four -point  form  of  these  formulas  appears  most  appropriate  for  the  present 
problem. 

First  consider  the  general  situation  where  the  solution  is  known  up  to 
qr  and  the  problem  is  to  find  the  values  of  9  and  y  and  their  derivatives 
at  rjr  +  ^  Tj  =  rjr+1  by  use  of  equations  (6.39  to  6.1-2)  and  equations  (6.49, 
6.51  and  6.53).  A  special  procedure  will  be  required  to  get  started  near  the 
wall,  and  it  is  described  in  succeeding  paragraphs.  Considering  first  the 
momentum  equation,  the  integration  indicated  in  equations  (6.39,  6.4-1,  6.42) 
will  be  approximated  by  the  extrapolation  and  interpolation  formulas.  First 
the  extrapolation  formulas  use  values  of  (C  cp")'  and  9"  at  the  r,  r-1, 
r-2,  r-3  stations  to  determine  values  of  cp,  cp',  9"  and  (C  9")'  at  the 
r+1  station.  The  formulas  are 

(C  *  >r+l)  =  '  )r  +  |  55(C  9' '  )'r  -  59(C  9'  ' 

+  37(c  9”  ^-2  -9(c  9,,)^_5  ]  (6.60) 

where  the  subscript  E  is  used  to  denote  extrapolation.  The  step  length 
is  a  constant  in  the  solution.  The  error  in  (6.60)  is 

\  H  U  n)5  (c  9"U)]V 
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The  extrapolation  formulas  for  cp'  and  CP  are: 


^D.  -*i*& (»  C  “  59  Ci  *  37  Ca  "  9  C3>  (6-6l) 


with  an  error  of 


\  S  +  (an)5  9^(6) 


and 


Vl)E  ■  *r  +  ^  »J  + 


2  r 


with  an  error  of 


323  <P”  -  264  <p^!  +  159  9”2  “  3s  %] 

(6.62) 

Ej  s  +  ^  (An)6  ,«(») 


The  value  of  (C  <p  *  * )  at  r+1  can  now  be  determined  by  using  the  momentum 
equation  (6.37)  and  the  extrapolated  values  of  CP***  cp',  and  cp.  It  will  be 
denoted  by 

<° *' '  WjE  -  VC>  **W 

The  interpolation  formulas  can  now  be  used  to  determine  more  exact  values  of 
cp",  cp',  CP  and  (C  cp" )  at  the  r+1  station.  The  formulas  are 

(C  cp'  •  )r+1  -(C  ?")r  +  |  [9(C  cp*  •  )^+l)E  +  19(C  <p"  )'r  -  5(C  CP*  *  )j_x 

+  (c9**)^2]  (6*63) 

with  an  error  of 

s-^(an)5  (c  *••(») ]v 

and 


cp'  «  cp*  + 
Yr+1  Tr 

with  an  error  of 

and 


+  19  <P” 


K'-l 


CP' 

r 


-2 


85  s  (An)5  cpVI(t) 


(6.64) 
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(6.65) 


Vl  =  ^r  +  *r 


-  $bO 


58  9”x  +  171  <P“  -  36  q>^L  +  7  <P”2 


with  an  error  or 

E6  S  “  lSo  CPVI(?) 

Finally  the  value  of  (C  cp' '  is  obtained  by  again  using  the  basic  momentum 

equation  (6.57)  and  the  interpolated  values  above 

(C  cp*  •)*  =  F1(cp",  cp«,  9) 

A  comparison  of  the  error  terms  for  the  extrapolated  and  interpolated  values 
shows  that  not  only  are  the  interpolation  errors  much  less  than  the  extra¬ 
polation,  they  are  also  opposite  in  sign.  Therefore  the  exact  value  of  the 
quantity  in  question,  say  cp',  must  lie  in  the  bounds  of  the  extrapolated 
and  interpolated  values .  Thus  there  is  a  check  on  the  accuracy  of  the  pro¬ 
cedure.  This  method  of  checking  is  the  concept  used  in  all  predictor- cor¬ 
rector  methods  of  numerical  analysis  (for  example,  Milne's  methods  described 
in  Reference  12,  pages  199-202).  The  method  developed  here,  losing  the  Falkner 
and  Adams  formulas,  was  chosen  over  the  more  established  predictor-corrector 
methods  because  its  use  of  higher  derivatives  gives  much  smaller  errors.  It 
is  particularly  well  suited  to  solving  high-order  ordinary  differential 
equations.  The -solution  can  be  made  as  exact  as  desired  by  choosing  a  small 
enough  step  size  Aq. 

The  use  of  both  the  extrapolation  and  interpolation  formulas  does 
essentially  mean  that  the  integration  over  each  step  length  is  being  per¬ 
formed  twice.  In  the  interest  of  saving  computing  time,  a  study  was  made  to 
see  if  use  of  only  the  extrapolation  formulas  would  give  accurate  solutions. 

In  some  flows  it  was  found  that  the  extrapolation  gave  oscillating  values  of 
9"  and  9'  at  the  outer  edge  of  the  boundary  layer.  These  oscillations 
disappeared  when  both  the  extrapolation  and  interpolation  formulas  were  used. 
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ftie  formulas  for  performing  the  Integrations  required  In  solution  of  the 
energy  equation  are  similar  to  the  shove.  They  are: 

for  equation  (6.51) 

VlJj.  "  Tr  +  $  [55  “  59  K-l  +  37  K-2  ~  9  ^-3  ]  (6*66) 

for  equation  (6.53) 

♦r.1), -*,4  [55  ♦;  -  59  ♦;.!  ♦  37  »;.2  -  9  ]  (6.67) 

and  from  equation  (6.49) 

Tr+1)B  “y2^V  *eW 
The  interpolated  values  are  given  by 

Vi  ■  'r  ♦  f-  [<9  *■»!),  *  W  K  -  5  ’’J-i  +  Va  ] 

♦r+1  -  *r  *  §  [  9  ♦»!  '»*;->  *;.l  +  ♦i.a  ] 

and  finally  from  equation  (6.4-9) 


ir'  .  =  F_(ir.  a)  . 


The  extrapolation-interpolation  formulas  above  require  values  of  the 
variables  at  four  previous  tj  stations.  To  get  started  at  the  wall  Taylor's 
series  will  be  used.  Consider  the  Taylor's  series  for  C  9" 

2 

(c  <P")An  =  (c  9")w  +  ATj(C  9,,)v  +  (C  <P")”  +  ... 

The  use  of  high-order  terms  in  the  expansion  would  require  values  of  the 
derivatives  of  C  at  the  wall,  for  example,  the  use  of  the  third  term  re¬ 
quires  the  value  of 


(6.68) 

(6.69) 
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1 1  Ty 

(C  q>")w  -  (C"<p")w  +  2(C'<p‘")w  +  (C  <p  )w 

Derivatives  of  C  are  not  knovn  from  the  method  of  solution,  and  it  would 
be  difficult  to  obtain  their  values  accurately.  Therefore  only  two  terms 
will  be  used  in  the  Taylor's  series  expansion.  Error  studies  show  that  if 
the  same  accuracy  is  to  be  maintained  in  the  expansion  as  in  the  extrapola¬ 
tion  —  interpolation  formulas  above,  much  shorter  steps  in  q  are  required 
in  the  Taylor's  series.  The  study  shows  that  for  about  five-place  accuracy 
in  the  values  of  cp  and  a  step  size  of  A?)  =*  0.1  is  sufficient  in 
the  extrapolation  —  interpolation  formulas  but  that  the  Taylor's  series 
expansion  requires  a  step  length  of  about  0.006.  The  accuracy  can  be  held 
by  using  the  Taylor's  series  to  obtain  the  values  at  only  T  =  ^rom 
the  wall.  Here  AH  is  used  to  denote  the  step  size  used  in  the  four-point 
extrapolation  —  interpolation  formulas  above.  Values  are  then  built  up  to 
the  full  length  step  AH  by  using  two-point  and  three-point  extrapolation 
formulas.  Accuracy  requires  a  step  size  of  at)/4  ?or  the  two-point  form, 
and  a  step  size  of  An/2  for  the  three-point  form.  The  steps  in  the  pro¬ 
cedure  are  outlined  in  Table  I  and  are  sketched  in  figure  6.  The  equations 
to  be  used  are: 


TAYLOR*  SHOES  -  STEP  SIZE  dfl/l6 


*  = 

(C  <p**)w  +  f|  (C  <P**)J 

(6.70) 

16  = 

-1+S 

(6.71) 

%/l6  = 

.  —  :*'l  +  V- 1  .11 

16  512  ¥W 

(6.72) 

W16  " 

\  +  S*s 

(6.73) 

W16  * 

*» +  H  K 

(6.7*) 
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TWO-POINT  FORMULAS  -  STEP  SIZE  &t/l6 


<c  ’"Wa  - 

(0  *"Wl6  +  f  [3(0  -  <c  *">w  ] 

(6-75) 

II 

1 

e- 

’’ww ♦  f  [3  »y/l6  -  »;•] 

(6.76 

%/8  " 

2 

%/l6  +  §  ’i/16  +  I5?T"  [4  ’i/16  “ 

(6.77) 

%/8s 

%/l6  +  ft  [3  *i/l6  ~  K  ] 

(6.78) 

♦/»l/8  = 

Wi6  +  f  [3  ♦An/ie  “  »v  ] 

(6.79) 

TWO-POINT  FORMULAS  —  STEP  SIZE  A»j/8 

<°  *MWi.  - 

(c  ’"Wa'f  [5(c  ’",)W8-(C  ] 

(6.80) 

9V>' 

’Vs + 1?  [5  »y/s  -  *;•  ] 

(6.8l) 

H 

\ 

/ 

2 

%/8  +  f1  ’i/8  +  [k  ’i/8  “ 

(6.8e) 

%A  = 

%/8  +  %  [5  Ti/8  “  K  ] 

(6.83) 

♦ah/4 

Ve  +  i?[5  ♦i/s  “  ♦*  ] 

(6.84) 
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TWO-POUT  FORMULAS  -  STEP  SIZE  tsjk 


<c  »*')rtl  ■  (C  ^  [~3<C  V")'r  -  (c  *••),.!  1  (6.05) 

♦«1  *  +  f1  [3  »"  -  %  ]  (6-86) 

Vi  ■  »r  ♦  f1 *; ♦  [*  *;■-  »"i  ]  <6-87) 

*r+l  ■  *r  *  f1  [5  *?  -  TLl  ]  (6-88) 

♦r+1  ■  ♦,  +  f  [5  K  -  K-l  ]  (8-89) 


THRKE-POUTT  FORMULAS  —  STEP  SIZE  Aij/2 


(C  *">r+l  =  <c  *">r  +  |J  [2MC  -  l6(C  f")^  +  5(C  9»*)^2  ] 

^r+l  =  +  §  [25  -  16  t*!!  +  5  <p*_'2  ] 

19  -  10  9*^  +  3  q»»*2j  (6.92) 


(6.90) 

(6.91) 


r+1 


r+1 


*r+l 


\  *  If  [«  K  -  16  K-l  +  5  r;_2  j 

♦r  +  f?  [2J  *;  -  16  ♦;.!  +  5  »;.a  ] 


(6.93) 

(6.96) 
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6.9  Outline  of  Procedure  for  Solving  Momentum  and  Baergy  Equations 

Simultaneously 

The  procedure  for  solving  the  momentum  and  energy  equations  Is  briefly  des¬ 
cribed  In  Section  6.k ,  and  the  reason  for  its  choice  is  discussed  there.  Details 
of  the  procedure  are  given  in  this  section.  Consider  the  general  case  when  the 
program  is  solving  the  equation  at  the  n-station.  Values  of  cp  and  f  and  their 
derivatives  at  all  previous  stations  will  be  known.  First  the  mcmentun  equation 
is  solved  using  the  fluid  properties  from  the  n-1  station.  The  values  of  q>'s 
from  this  solution  are  used  to  solve  the  energy  equations  but  still  using  the 
fluid  properties  from  the  n-1  station.  Then  new  fluid  properties  are  determined 
and  an  iterative  procedure  followed  until  convergence  of  solutions  of  both  the 
momentum  and  energy  equations  is  obtained.  In  the  iterative  procedure  let  q  =  i 
(an  integer)  indicate  one  solution  of  the  momentum  equation  with  the  accompanying 
solutions  for  enthalpy  and  fluid  properties.  The  procedure  is: 

SLl. 2 

(a)  The  momentum  equation  is  solved  by  using  the  fluid  properties  from 
the  n-1  station.  It  is  solved  by  the  cut-and-try  and  interpolating 
procedure  described  in  Section  6.5.  The  solution  is  denoted  as  cpQ. 

(b)  The  energy  equation  is  solved  by  using  the  cpQ  value  and  the  n-1 
fluid  properties.  The  solution  is  denoted  as  (♦o)jj=0/  where  L  is 
a  count  of  successive  solutions  of  the  energy  equation  and  fluid 
properties  at  a  given  q. 

(c)  The  solution  (*  )T  „  is  used  to  determine  new  fluid  properties. 

They  sure  denoted  as  (FP^I^O' 

(d)  The  fluid  properties  (FP  )x__  and  the  solution  ©  are  vised  to 

O  o 

obtain  a  second  solution  to  the  energy  equation  and  new 

fluid  properties  (FP  )T 

O  JjsJ. 

(e)  This  iterative  procedure  is  continued  until  L  *  L-MAX.  Proceed 
to  q  =  l. 

3jlA 

(a)  The  momentum  equation  is  solved  a  second  time  by  using  fluid  proper¬ 
ties  ( FPq  )l_max  8214  "the  solution  is  denoted  as  <p^. 

(b)  Steps  (b)  through  (e)  in  q  =  0  are  repeated  to  obtain 
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9l>± 

(a)  The  procedure  in  Q  =  1  is  repeated,  using  always  the  latest  values 
of  (9)#  (♦),  and  (FP)  until  either 

Q  =  SttX 

or 

‘CIq-'CVi  <  * 

where  c  is  an  accuracy  input.  If  either  of  these  conditions 
is  satisfied,  the  program  proceeds  to  the  next  station. 

The  procedure  is  sketched  in  figure  7,  and  further  details  of  the 
computer  program  are  given  in  Appendix  B. 

6.10  Starting  the  Solution 


At  x  =  0  the  x-dependent  terms  disappear  in  both  the  momentum  and 
energy  equations.  The  equations  reduce  to: 


3? 


(c 


9**)  =  -  T  —  C  <?'•  +  C  pf cp*2  +  2  9'  +  1 - — 1  -  C  H(cp  +  Tj)cp* » 

P  00  [_  P  J  00 


(6.95) 


for  momentum  and 


6 

3n 


^  C(1  )(cp.  +  i)q>» 

e 


-  T 


£  r  +  -jp  c(l  -  ^  )(9'  +  1)  9'  *J  -  cmn(9  +  n)  ♦'  (6.96) 


for  energy.  Hence  values  of  9  and  at  previous  stations  are  not  required. 
The  x-dependent  terms  also  disappear  for  similar  flows.  Such  flows  occur  when 
P,  R  and  the  wall  boundary  conditions  are  constants  for  all  x's.  They  In¬ 
clude  those  with  constant  pressure,  such  as  flows  about  a  flat  plate  and 
wedges.  Similar  -flows  are  discussed  more  completely  in  Chapter  8  of  Hayes 
and  Probstein,  Reference  13. 
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Figure  7  -  Flow  diagram  for  solving  boundary-layer  equations  at  station  n. 
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The  procedure  of  solution  requires  that  the  momentum  equation  be  solved 
first,  but  to  do  so  requires  values  of  the  fluid  properties.  At  a  downstream 
station  the  properties  are  approximated  for  the  first  solution  of  the  momentum 
equation  by  assuming  them  to  be  the  same  as  those  at  the  previous  station.  To 
get  started  at  the  x  =  0  station,  a  linear  enthalpy  profile  that  satisfies 
the  inner  and  outer  boundary  conditions  for  enthalpy  is  assumed.  The  fluid 
properties  obtained  from  this  enthalpy  profile  are  then  used  to  start  the 
solution.  After  the  first  solution  of  the  momentum  equation  at  any  station 
is  found,  the  fluid  properties  obtained  from  the  latest  solution  of  the  energy 
equation  are  used  for  subsequent  solutions  of  the  momentum  equation.  The 
iterative  procedure  is  presented  in  the  preceding  section  of  this  report. 


Values  of  P  and  R 

The  values  of  P  and  R  are  either  inputs  in  the  program  or  are  cal¬ 
culated  from  the  velocity  distribution  ug  versus  x,  and  the  radius  distri¬ 
bution  rQ  versus  x.  They  are  defined  by 


(6.13) 


x  drp 
rp  dx 


1  +  T 


dn 


(6.1*0 


Both  Fn  and  (rQ/x)(dro/dx)  are  always  inputs  at  x  =  0,  but  may  be  calculated 
at  aft  stations  by  use  of  Iagrangian  derivative  formulae 


and 


=  An 


.  u  +  A  u  +  A  ,  u 
-1  en-l  n  en  n+1  en+l 


A  i  r  +  A  r  +  A.  r 
n-1  °n-l  n  °n  n+1  °n+l 


(6.97) 


(6.98) 


where 


Vl  "  lJShl-V(*aJ.“W 


2  x  —  x  ,  —  x  , 

.  n  n+1  n-1 

n  ' 

.  _  n  n-1 _ 

«+l  "  (x  ,,  —  X  ,  )(X  ,,  —  X  ) 
'  n+1  n-1  '  n+1  n' 


(6.99) 


6.11  Boundary-Layer  Parameters 


Once  the  profiles  of  <p  and  \|r  and  their  derivatives  have  been  de¬ 
termined  at  an  x-station,  the  program  determines  the  conventional  boundary- 
layer  parameters  of  displacement  thickness,  momentum  thickness,  local  skin 
friction  and  local  heat  transfer.  First  the  profiles  of  cp,  and  their 
derivatives  are  transformed  to  the  more  conventional  profiles  of  the  stream 
function  f,  the  velocity  ratio  f,  enthalpy  ratio  g,  etc. 


The  displacement  thickness  is  calculated  from  the  well  known  formula 


# 

& 


(6.100) 


In  the  true  definition  of  a  displacement  thickness , equation  (6.100)  is  exact 
only  for  two-dimensional  flow.  For  axially  symmetric  flow  the  exact  dis¬ 
placement  thickness  is  given  by  the  quadratic  equation 


8  (1 


* 

6 

2r~~ 

o 


00 

cos  a)  =  Ml  +  -X—  cos  a)(l  —  ft  )  dy 
v  ro  pee 


(6.101) 


Because  of  difficulties  in  solving  (6.101),  the  simpler  expression  (6.100)  is 
used  to  calculate  displacement  thickness  for  all  flows.  Transforming  to  the 
x,  7)  coordinates,  (6.100)  is 
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(6.102) 


The  program  calculates  a  dimensionless  displacement  thickness  A  ,  which  is 
defined  by 


u 


08  X 


u  c 
e 


(6.103) 


When  comparing  solutions  with  real  flows  the  dimension  tj  is  not  con¬ 
venient,  because  of  the  transformation  between  it  and  the  true  physical  height 
y.  Therefore  the  program  calculates  y  from  tj. 


(6.4) 


Define  a  dimensionless  Y  as 


Then 


(6.104) 


(6.105) 


(6.106) 
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Introduction  of  Y  in  (6.103)  gives 


* 

A  =  lim 
H— eoo 


L  V  ue  c  Pe  J 


(6.107) 


The  momentum  thickness  is  given  by 


(6.108) 


or  in  x,  tj  coordinates  by 


0 


w 

J f*(l  -  f') 


dl 


(6.109) 


The  program  calculates  a  dimensionless  momentum  thickness  defined  as 


„  Ip  UT 

0-^  I 


(6.U0) 


or 


(6.111) 


The  shear  stress  at  the  wall  is 


/6u\ 

Tw=  *Sr  (3yV 


(6.112) 


or  in  x,  t)  coordinates 


/  u" 

/ _ e 

Tw  "  \/p  u 
V  'ce  oc 


e  w 


(6.113) 


Introduction  of  the  local  skin- friction  coefficient  defined  as 
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(6.114) 


cf  aT 


T 

2  PCBU«8 


gives 


=  2^  (Ji) 

V  PcoU«X  \  Uco  / 


3/2 


pw 

p»  ^ 


w 


A  conventional  shear  parameter  is 


* 

c  =  c 
cf  f 


vs 


c  /  u  i 

*  l\ 


\  ^ 


=  2 


UJ 


PW 


ft. 


The  heat  transfer  at  the  wall  is 


n  -  /  k  6h 

^ =  1  °p  ^  , 


Now 


6h 


_  3 

“~57 


3H1 


w 


since  =  0.  Heat  transfer  in  x,  T]  coordinates  is  given  by 


^  VpJV1 


“w  He 


The  Stanton  number,  a  heat  transfer  parameter,  is  defined  as 


St 


_  ~~  Sf  _  /  ^oo  /  ue  Pw^W  1  ®w 

=  (Hg  -  Hw)  “  V  p^ux  \J  U{o  p^  Prw  1  -  ^ 


The  ratio  of  Stanton  number  to  skin  friction  is  given  by 


St 


2  Pr 


u-  i  < 

ue 


*  The  Stanton  number  is  most  often  written  in  terms  of  H _ —  H 

rec  w 

than  Hg  —  Hw,  but  this  is  not  convenient  in  the  present  method. 


(6.115) 

(6.116) 

(6.117) 

(6.118) 

(6.119) 

(6.120) 

(6.121) 

rather 
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7.0  CHARACTER  OF  SOLUTION  AS  LEARNED  FROM  TRIAL  RUNS 


The  principal  purpose  of  this  section  is  to  establish  the  accuracy  and 
character  of  the  method  of  solution.  It  will  be  done  by  presenting  a  large 
number  of  flow  cases  that  include  large  ranges  in  pressure  gradient,  heat 
transfer,  transverse  curvature  effects,  and  different  fluid-property  laws. 
Whenever  possible,  comparisons  are  made  with  exact  solutions  or  with  experi¬ 
ment.  Also,  in  some  cases  comparisons  are  made  with  approximate  methods  of 
other  investigators.  Section  7.1  discusses  the  accuracy  for  similar  flows. 

Since  these  are  almost  the  only  flows  for  which  exact  solutions  are  available, 
they  are  used  to  study  the  effect  of  the  various  inputs  in  the  computer  pro- 
•gram  on  the  accuracy  of  the  calculations  in  Section  7.2.  Section  7.3  compares 
the  present  method  of  solution  with  that  of  Chapman  and  Rubesin  for  the  special 
case  of  a  flat  plate  with  variable  wall  temperature.  Effect  of  the  various 
fluid-property  laws  on  recovery  factor  is  determined  in  Section  7.^*  Section 
7«5  compares  the  present  method  with  the  finite-difference  method  of  Fltfgge- 
Lotz  and  Blottner.  The  present  method  is  compared  with  experimental  results 
for  the  special  case  of  a  circular  cylinder  in  Section  7*6.  An  example  of 
internal  flow,  the  boundary  layer  inside  a  nozzle,  is  presented  in  Section  7.7. 
Boundary  layers  on  a  reentry- type  body  sire  presented  in  7.8  showing  the  effect 
of  altitude ,  body  temperature,  and  the  fluid-property  laws  on  the  solutions. 
Transverse  curvature  effects  are  studied  in  7*9.  Finally,  in  Section  7.10, 
discontinuities  in  body  temperature  and  mass  transfer  are  studied.  This  study 
is  applicable  to  some  types  of  ablation  cooling. 

The  first  —  and  an  obvious  —  check  of  the  program  was  simple  incompressible 
flow.  For  such  flow  the  equations  reduce  to  the  form  presented  in  References 
1  and  2,  where  the  authors  developed  the  present  method  of  solution  for  in¬ 
compressible  flow.  Many  of  the  cases  presented  in  References  1  and  2  were 
recalculated  by  the  present  method.  Results  were  identical  to  the  earlier 
ones  and  will  not  be  repeated  here.  Since  the  energy  equation  is  negligible 
and  fluid  properties  are  constants  for  incompressible  flow,  these  cases  afforded 
a  check  on  only  -that  portion  of  the  program  that  solves  the  momentum  equation. 
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7.1  Similar  Flow 


Since  almost  the  only  flows  for  which  highly  accurate  solution  are  known 
are  the  similar  flows  calculated  by  Cohen  and  Reshotko  (Reference  14),  these 
were  first  used  to  check  the  programing  of  the  equations.  Similar  flow  occurs 
when  the  solution  is  a  function  of  tj  only,  that  is,  when  f  and  g  and  their 
derivatives  are  independent  of  x.  In  their  method  Cohen  and  Reshotko  used 
Stewartson's  transformation  to  transform  the  compressible  equation  to  an  in¬ 
compressible  form.  Specifically,  the  transformations  are: 

x 


P- 


dx 


(7.1) 


(7.2) 


(7.3) 


where  here  X,  Y,  and  U  are  the  transformed  quantities,  the  subscript  T 
refers  to  free-stream  stagnation  values,  a  is  the  velocity  of  sound,  and 
Tv  is  a  constant  in  a  linear  viscosity- temperature  relation.  The  transformations 
are  restricted  to  the  cane  where  viscosity  varies  linearly  with  temperature  and 
for  Prandtl  number  1.0.  The  same  assumptions  are  used  in  the  present  method  when 
comparing  its  results  with  those  of  Cohen  and  Reshotko.  Similar  flow  occurs 
when 


U 


(7.4) 


Cohen  and  Reshotko  present  their  results  for  specific  values  of  the  pressure 
gradient  parameter  fi,  which  is  related  to  m  by 


P  = 


2m 

m  +  1 


(7.5) 


Results  of  the  present  method  are  compared  with  those  of  Cohen  and  Reshotko 
in  Table  II  for  various  values  of  m  and  g^.  In  the  present  method  the  calcu¬ 
lation  can  be  compared  directly  with  Cohen  and  Reshotko's  results  by  letting 


(7.6) 


u  =  U  =c.  x” 
eel 

That  is,  essentially  low  speed  flow  is  assumed,  but  density  does  vary  with 
enthalpy 

pe  h 

T  =  “"K“ 

The  parameter  P  is  identically  equal  to  m  for  such  flows.  Results  of  the 
two  methods  differ  by  only  2  in  the  fourth  decimal  place  for  favorable 
pressure  gradients.  Cohen  and  Reshotko  presented  four  significant  figures  in 
their  results.  In  the  present  method  the  various  inputs  that  affect  the  accuracy 
of  the  results  were  selected  so  that  the  calculations  converged  to  five  signifi¬ 
cant  figures.  These  inputs  are  discussed  in  the  next  section.  The  agreement 
between  the  two  methods  is  not  so  good  near  separation,  that  is,  near  f ' 1  =0. 

In  the  worst  case,  the  difference  is  2  in  the  second  decimal  place  in  fj'. 

The  present  method  becomes  very  sensitive  near  separation  (fj 1  <  0.08)  and  its 
difficulties  in  obtaining  a  solution  here  are  discussed  in  detail  in  References 
2  and  3  for  the  case  of  incompressible  flow.  The  values  of  and  fj ' 
presented  in  Table  II  for  P  =  20  and  P  =  2.0  should  not  agree  exactly.  A 
value  of  p  =  2.0  corresponds  to  a  value  of  P  =  oo ,  and  it  is  of  course  im¬ 
possible  to  calculate  such  a  flow  by  the  present  method.  A  value  of  P  =  20 
was  used  to  convert  the  values  of  f^'  and  presented  In  Reference  1^  for 
P  =  2.0  to  the  values  presented  in  Table  II.  As  P  becomes  very  large  the 
results  of  the  present  method  approach  those  of  Cohen  and  Reshotko  for  P  =  2.0. 


For  the  case  of  favorable  pressure  gradients  with  heated  walls,  that  is, 
the  wall  temperature  greater  than  recovery  temperature,  the  velocity  within 
a  portion  of  the  boundary  layer  can  exceed  the  local  external  velocity.  This  is 
known  as  velocity  overshoot.  An  example  is  shown  in  figure  8,  where  P  =  l/3, 
gw  =  hw/Hg  =  2.0,  and  =  3.0.  Both  enthalpy  and  velocity  profiles  are  shown. 
The  assumptions  for  fluid  properties  are  the  same  as  those  used  by  Cohen  and 
Reshotko;  Pr  =  1.0  ana  a  viscosity  proportional  to  temperature.  By  means  of 
the  Stewartson  transformation  the  flow  reduces  exactly  to  that  presented  in 
Reference  14  for  the  transformed  incompressible  flow  with  P  =  l/3  and  g  =  2.0. 
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TABLE  II 

COMPARISON  OF  SHEAR  PARAMETER  AND  HEAT-TRANSFER  PARAMETER  CALCULATED 
BY  PRESENT  METHOD  WITH  THOSE  CALCULATED  BY  COHEN  AND  RESHOTKD 


f " 
w 

K \ 

0 

P 

Sv 

Cohen  and 

Present 

Cohen  and 

Present 

Reshotko 

Method 

Reshotko 

Method 

1.0 

1.0 

1.0 

1.0 

1.0 

2.0 

1.23259 

1.73671 

t 

0 

o\ 

H  O 
£ 

0 

-0.61533 

0.50 

0.33333 

0.0 

0.4741 

0.47413 

0.4040 

0.40399 

0.50 

0.33333 

0.2 

0.5346 

0.53477 

0.3290 

0.32922 

0.50 

0.33333 

0.6 

0.6489 

0.64907 

0.1706 

0.17068 

0.50 

0.33333 

1.0 

0.7575 

0.75745 

0 

0 

0.50 

0.33333 

2.0 

I.OO85 

1.00863 

-0.4674 

-0.46745 

0.0 

0.0 

0.0 

0.3321 

0.332057 

0.3321 

0.332057 

0.0 

0.0 

0.2 

0.3321 

0.332057 

0.2657 

O.265676 

0.0 

0.0 

0.6 

0.3321 

0.332057 

0.1328 

0.132838 

0.0 

0.0 

1.0 

0.3321 

0.332057 

0.0 

0.0 

0.0 

0.0 

2.0 

0.3321 

0.332057 

-0.3321 

-0.332057 

-0.10 

-0.04762 

mu 

0.1246 

0.12478 

-0.2783 

-0.27826 

-0.160 

-0.07407 

Wwl 

0.1296 

0.12809 

0.0 

0.0 

-0.1988 

-0.09041 

1 

0.0 

0.02405 

0.0 

0.0 

-0.200 

-0.09091 

0.1472 

0.14747 

0.1096 

0.10873 

-0.240 

-0.10714 

0.6 

0.07H 

0.08172 

0.0985 

0.10304 

-0.300 

-0.13044 

0.0 

0.2098 

0.21213 

0.2810 

0.28392 

-0.300 

-0.13044 

0.2 

0.1376 

0.13843 

0.2080 

0.21079 

-0.325 

-0.13979 

0.2 

0.0888 

0.09720 

0.1910 

0.19981 

* 

2.0 

* 

20.0 

0.0 

2.3917 

2.33313 

1.6860 

1.62839 

2.0 

20.0 

0.2 

2.0719 

3.01625 

1.4034 

1.38057 

2.0 

20.0 

0.6 

4.3191 

4.23688 

0.7466 

0.74326 

2.0 

20.0 

1.0 

5 . 4666 

5. 3^950 

0 

0 

2.0 

20.0 

m 

8.0615 

7.89056 

-2.1429 

-2.12336 

Since  6  =  2.0  corresponds  to  P  =  od  ,  the  two  methods  can  not  agree 
exactly . 


7.2  Effect  of  the  Various  Inputs  on  the  Accuracy  of  the  Solution 


Since  exact  answers  are  available  for  similar  flows,  these  flows  have 
been  used  to  study  the  accuracy  of  the  computer  program.  Quantities  that 
are  inputs  and  that  affect  the  accuracy  of  the  computer  results  are  A^, 

K,  L,  e,  AX,  and  the  trial  values  of  either  gj  or  g^. 

The  step  length  in  rj  affects  the  accuracy  of  the  calculation  because 

of  truncation  errors  in  the  integration  formulas.  Effect  of  the  step  length 

on  the  calculation  of  f^'  is  shown  in  Table  HI.  Its  effect  on  gj  is 

similar.  In  fact  the  error  in  f"  is  typical  of  the  maximum  error  that 

w 

occurs  in  all  quantities  calculated,  including  the  enthalpy  profiles.  In  the 
table  and  in  the  results  that  follow,  the  various  inputs  that  effect  the 
accuracy  have  the  following  values ,  unless  otherwise  specified:  ah  =  0.05, 

=  9*0,  K  =  0.50,  L  =  3.0,  e  =  0.000001.  The  particular  values  of  P 
and  gtf  presented  in  the  table  were  chosen  because  exact  results  (4  decimal 
places)  are  available  for  them  (Reference  14)  and  because  they  cover  a  large 
range  of  both  pressure  gradient  and  wall  temperature .  The  table  shows  that 
for  favorable  pressure  gradients  a  At)  of  0.05  gives  five-place  accuracy  and 
a  Arj  of  0.2,  three-place;  but  the  error  is  larger  for  unfavorable  pressure 
gradients  (P  <  0). 

The  maximum  value  of  rj,  that  is,  t^,  is  an  input  in  the  program.  The 
method  of  solution  (see  Sections  6.5  and  6.6)  forces  the  boundary  layer  to 
satisfy  the  external  conditions  on  <p'  and  \|r  at  this  height.  Boundary- 
layer  thickness  can  vary  greatly  with  Mach  number,  pressure  gradient,  and  wall 
temperature,  gy.  If  is  chosen  too  small,  the  results  will  be  in  error. 

Effect  of  different  values  of  on  two  of  the  similar  flows  is  shown  in 

Table  IV.  Values  of  f ' ’  are  shown,  but  the  effect  on  g^  is  similar. 

Because  can  vary  greatly  with  the  flow  conditions,  especially  at  high 
Mach  numbers ,  it  is  not  always  possible  to  select  its  proper  value  without 
first  making  a  trial  run.  It  must  be  selected  so  that  the  higher  derivatives 
of  f '  and  g  approach  zero  asymptotically  at  the  outer  edge  of  the  boundary 
layer.  In  all  the  cases  that  have  been  calculated,  it  was  found  that  the 
maximum  value  of  for  five-place  accuracy  was  about  9*  For  favorable 
pressure  gradient,  can  be  considerably  less  for  the  same  accuracy. 
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TABLE  III 

EFFECT  OF  STEP  SIZE  IN  T)  ON  ACCURACY  OF  SOLUTION 
=  9*0,  K  =  0.5,  L  =  3,  e  =  10”^ 


FLOW 

VALUES  OF  f'« 
w 

P 

Sv 

£ 

II 

O 

• 

\Ji 

An  =  0.2 

ESI 

ah  =  0.05 

An  =  0.02 

1.0 

0.333333 

-0.0476191 

-0.1304348 

1.0 

0.2 

2.0 

0.2 

1.270* 

0.540 

0.1905 

_ ** 

1.23270 

0.53440 

0.1234 

O.169 

1.232627 

0.53690 

0.12391 

0.1391 

1.232591 

0.534770 

0.124780 

0.138425 

1.232588 

0.534786 

0.124702 

0.138256 

* 

Digits  underlined  may  not  be  significant 

Program  unable  to  converge  to  a  solution 

TABLE  IV 


EFFECT  OF  INPUT  n*  ON  ACCURACY  OF  SOLUTION 
AT)  =  0.05,  K  =  0.5,  L  =  3,  €  =  10"6 


f 

1 

P  =  1.0 

gw  =  1.0 

p  =  -0.0476191 
=  2‘° 

12 

1.232588 

0.124778 

10 

1.232588 

0.124780 

8 

1.232588 

0.12482 

6 

I.23259I 

0.1269 

4 

1.2326 

0.1342 

Digits  underlined  may  not  be  significant 

Selection  of  much  larger  than  necessary  can  affect  the  accuracy  of  the 
results,  for,  in  the  trial  solutions  of  the  energy  equation,  the  values  of 
t  and  T|r*  at  the  outer  edge  of  the  boundary  layer  can  grow  exponentially. 

In  cases  when  was  set  too  large,  values  as  high  as  10^  have  been  ob¬ 
served.  In  the  procedure  of  solution,  the  two  trial  solutions  are  added  to 
meet  the  outer  boundary  condition,  ^(r^)  =0  (see  equations  6.56  and  6.57). 
This  process  of  addition  can  lead  to  large  roundoff  errors  in  the  "correct" 
solution  of  ijr  if  the  trial  values  of  t  are  very  large.  (The  computer 
program  carries  only  eight  significant  digits.) 

The  bound  K  also  affects  the  accuracy  of  the  results.  Its  purpose  is 
described  in  detail  in  Section  6.5;  briefly,  it  affects  the  solution  of  the 
momentum  equation  in  the  following  way:  the  correct  value  of  f"  is  found 
by  trying  different  values  of  f  ^ ' .  Once  three  trial  solutions  have  been 
found  such  that  9'  at  is  between  the  bounds  of  —  K  S  cp 1  ( )  g  K,  a 

three-point  interpolation  procedure  is  then  used  to  determine  the  solution 
that  satisfies  the  outer  boundary  condition.  Effect  of  the  bound  on  the 
solution  is  shown  in  Table  V  for  two  similar  flows. 


TABLE  V 

EFFECT  OF  THE  BOUND  K  ON  ACCURACY  OF  SOLUTION 

-6 

ATI  =  0*05;  T*  *  9.0,  L  =  5,  6  =  10 


K 

f  " 
w 

P  =  1.0 

gw  =1.0 

p  =  -0. 01+76191 

v  2-° 

0.02 

1.232588 

0.121+780 

0.2 

1.232588 

0.121+780 

0.5 

1.232591 

0.121+786 

1.0 

1.23260 

0.121+81 

2.0 

1.23262 

0.12l+8£ 

Digits  underlined  may  not  be  significant 
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Results  of  a  large  number  of  cases  show  that  a  K  of  0.5  produces  values 
accurate  to  five  decimal  places,  providing  the  Mach  number  is  not  so  high  and 
the  wall  so  warm  as  to  give  large  overshoot.  In  such  extreme  cases  the 
velocity  ratio  f '  =  1  +  cp'  can  exceed  1.5  in  the  boundary  layer.  In  the 
example  shown  in  figure  8,  f '  has  a  maximum  value  of  1.49.  In  such  cases 
K  must  be  chosen  larger  than  the  amount  of  overshoot.  But  in  favorable 
pressure  gradients  care  also  must  be  taken  that  K  is  not  too  large,  no 
more  than  about  2.0.  For  it  is  possible  in  such  cases  that  though  9*  for 
low  trial  solutions  at  first  grows  negatively,  the  growth  can  reverse;  cp' 
then  grows  positively  for  large  t),  and  the  low  trial  solution  can  exceed 
a  high  one  (see  Reference  3 ) •  If  this  happens,  the  program  fails  to  converge 
to  an  answer. 

The  input  L  determines  the  number  of  successive  solutions  of  the 
energy  equation  used  to  calculate  fluid  properties  and  the  enthalpy  profile 
(see  Section  6.9).  Remember  that  in  solving  the  energy  equation,  the  fluid 
properties  are  first  assumed  and  the  enthalpy  distribution  is  computed. 

New  fluid  properties  are  determined  from  the  calculated  enthalpy  and  the 
energy  equation  is  solved  again.  This  iterative  procedure  is  continued  L 
times.  From  the  cases  studied  it  has  been  found  that  for  flows  where  the 
enthalpy  changes  by  a  factor  of  about  10  across  the  boundary  layer  an  L  of 
3  produces  values  of  both  fluid  properties  and  enthalpy  that  converge  to 
at  least  four  decimal  places.  An  L  of  4  gives  convergence  to  at  least 
six  decimal  places.  When  the  enthalpy  change  through  the  boundary  layer  is 
much  less,  less  than  three  times,  an  L  of  2  gives  convergence  of  the 
calculated  values  to  at  least  four  decimal  places.  The  value  of  L  neces¬ 
sary  to  give  convergence  of  the  calculated  values  to  a  specified  accuracy 
depends  mostly  on  the  magnitude  of  the  enthalpy  change  across  the  boundary 
layer  and  is  independent  of  pressure  gradient,  Mach  number,  etc. 

In  the  method  of  solution,  once  the  energy  equation  has  been  solved 
corrected  fluid  properties  are  used  to  solve  again  the  momentum  equation  (see 
Section  6.9).  This  iterative  procedure  is  continued  until 

m'>a—  Vi|  <  •  ‘7.7) 
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where  c  is  on  input  and  Q  is  a  count  of  successive  calculations  of  the 
momentum  equation.  Again,  as  with  L,  the  number  of  times,  4  +  1,  that  the 
momentum  equation  has  to  be  solved  in  order  to  satisfy  eq. (7»7),  depends 
mostly  on  the  enthalpy  change  across  the  boundary  layer  and  appears  to  be 
independent  of  pressure  gradient,  Mach  number,  etc.  Variation  of  Q  with  c 
is  shown  in  Table  VI. 


TABLE  VI 

NUMBER  OF  SUCCESSIVE  CALCULATIONS  (Q  +  l)  OF  MOMENTUM  EQUATION 
NECESSARY  TO  SATISFY  ACCURACY  REQUIREMENT 


Enthalpy  Change 

Across 

Boundary  Layer 

<4 

c  =  10-6 

LT\ 

1 

O 

H 

II 

<* 

e  =  10~k 

€  =  10"5 

Approximately 

10  Times 

5 

3 

3 

2 

3  Times  or  Less 

3 

2 

1 

1 

The  effect  of  step  length  in  x  on  the  solution  appears  to  *be  identical 
to  that  reported  by  the  authors  in  References  2  and  3  for  the  incompressible 
boundary  layer.  The  details  will  not  be  repeated  here;  but,  briefly,  not 
only  must  the  size  of  Ax  be  considered  in  selecting  the  spacing  in  the 
x-directlon,  but  also  the  quantity  x/ax.  The  latter  appears  in  the  terms 
containing  the  x-derivatives  in  the  following  way.  From  the  momentum  equation 
the  term  containing  the  x-derivative  of  qp  *  is  replaced  by  a  finite-difference 
term  of  the  form 


x  <p 

* 


xn-xn-l/J 


x 

A  x 


cp*  a  <P* 
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As  ax  approaches  zero  AP1  approaches  zero  In  an  exact  solution;  but  in 
the  computer  program  roundoff  errors  exist  and,  furthermore,  these  errors  are 
multiplied  by  the  quantity  x/a  x.  Therefore  as  Ax  approaches  zero  the 
error  in  the  solution  is  not  necessarily  decreased.  This  type  of  roundoff 
error  is  observed  as  Ax  becomes  small  in  the  present  method.  It  does  not 
occur  for  values  of  x/a  x  less  than  about  25.  The  maximum  error  that  has 
been  observed  is  one  in  the  fifth  significant  digit  for  x/a  x  —  25.  The 
error  grows  as  x/a  x  becomes  larger.  Of  course,  this  roundoff  error  is  a 
function  of  the  number  of  significant  digits,  eight,  that  are  carried  in  the 
computer  program  and  is  not  inherent  in  the  method  of  solution.  It  can  be 
reduced  by  using  1 'double  precision* '  routines,  that  is,  sixteen  significant 
digits  would  be  carried  instead  of  the  present  eight.  In  all  of  the  calcu¬ 
lations  presented  in  this  report  A  x  was  selected  so  that  the  quantity 
x/a  x  was  less  than  25.  Therefore  it  is  believed  that  the  results  are  free 
of  this  type  of  roundoff  error. 

In  particular  cases  care  in  selecting  two  of  the  other  input  quantities 
is  necessary  to  prevent  erroneous  results.  If  gw  is  known, two  trial  values 
of  must  be  input;  or  if  g*  is  known,  then  two  trial  values  of  gw 
must  be  input.  In  the  method  of  solution  these  too  trial  values  are  used  to 
obtain  two  solutions  of  the  energy  equation  (see  Section  6.6),  Then,  because 
the  equation  is  linear,  the  correct  solution  is  found  by  combining  the  two 
solutions  to  meet  the  outer  boundary  conditions.  If  the  trial  values  at  the 
wall  are  far  from  the  correct  solution,  the  trial  values  of  y  and  i|r'  can 
become  very  large  near  the  outer  edge  of  the  boundary  layer.  This  is  the 
same  phenomenon  that  occurs  when  is  selected  too  large,  and  that  is  dis¬ 

cussed  above.  In  one  example  when  solving  for  which  had  a  correct  value 
of  approximately  0.4,  trial  values  of  0.0  and  1.0  were  input.  At  the 

the  outer  edge  of  the  boundary  layer  they  led  to  values  of  i|r  of  around 
9  10 

-10  and  +10  ,  respectively.  When  the  two  solutions  were  combined,  in  order 

to  satisfy  the  outer  boundary  condition  (see  equations  6.56  and  6 . 57 ) >  the 
large  trial  values  led  to  roundoff  errors  in  the  second  significant  digit 
of  i/.  This  roundoff  error  can  be  eliminated  by  selecting  trial  values  of 
g^  (or  gy)  near  to  the  correct  value,  within  one  significant  digit. 
Selection  is  done  by  either  experience  or  a  trial  run.  The  magnitude  of  the 
roundoff  error  can  be  easily  seen  by  inspecting  values  of  t  and  t|r*  near 
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7*3  Flow  on  a  Flat  Plate  with  Variable  Temperature 

Chapman  and  Rubes  in  developed  a  method  for  calculating  exactly  the  com¬ 
pressible  laminar  boundary  layer  on  a  flat  plate  with  variable  vail  tempera¬ 
ture  (Reference  15 ).  The  method  requires  that  the  surface  temperature  distri¬ 
bution  be  expressible  as  a  polynomial  in  surface  distance*  that  the  viscosity 
varies  linearly  with  temperature  through  the  boundary  layer*  and  that  the 
Prandtl  number  be  constant.  One  of  the  examples  presented  by  Chapman  and 
Rubesin  was  calculated  by  the  present  method  with  the  same  assumptions  for 
fluid  properties.  The  distribution  of  surface  temperature  is  given  by 

(0.25  -  0.83  X  +  0.33  X2)  (7.8) 

e,  e  e 

where  T^/T  =  1  +  0.169  is  the  wall  recovery  temperature  and  Tg  is 
the  temperature  at  the  edge  of  the  boundary  layer.  Since  constant  c^  is 
assumed  the  wall  enthalpy  distribution  has  the  seme  form  as  the  temperature. 
The  external  flow  has  a  Mach  number  of  3.0.  Results  of  the  two  calculations 
are  compared  in  figure  9  in  the  form  of 


versus  x.  The  latter  quantity  is  the  parameter  vised  by  Chapman  and  Rubesin. 
Greatest  disagreement  between  the  two  methods  is  5  in  the  fourth  significant 
digit.  Furthermore,  this  difference  is  consistent  for  the  whole  length  of 
flow;  the  values  calculated  by  the  present  method  are  (0.4  +  O.l)  percent 
less  than  the  values  calculated  by  Chapman  and  Rubesin.  In  the  present 
method  inputs  were  selected  in  such  a  way  that  the  calculation  was  accurate 
to  five  decimal  places.  In  using  the  method  of  Chapman  and  Rubesin  the  cal¬ 
culation  depended  on  parameters  given  to  four  figures  in  Reference  15. 

Therefore  the  two  calculations  might  be  expected  to  agree  to  at  least  four 
decimal  places .  But  Chapman  and  Rubesin  do  not  specifically  state  that  all 
four  digits  presented  are  significant.  The  value  of  recovery  factor  they 
calculated  is  presented  to  only  three  figures*  a  recovery  factor  <*  0.845. 

This  affects  the  values  of  the  other  quantities  —  temperature  and  heat  transfer— 
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calculated  by  their  method.  Recovery  factor  calculated  by  the  present  method 
for  the  flow  is  0.8475.  The  ratio  of  recovery  temperature  calculated  by  the 
present  method  to  that  by  Chapman  and  Rubes  in  is  1.0018. 

7.4  Effect  of  Fluid-Property  Laws  on  Recovery  Factor 


Effect  of  the  various  fluid-property  laws  on  recovery  factor  was  studied 
for  flow  over  a  flat  plate  and  flow  over  a  cone  at  Mach  numbers  from  1.0  to 
10.0.  Recovery  factor  rf  had  its  usual  definition 


1 


+  r. 


y  - 1 
2 


4 


(7.10) 


Relation  (7.10)  is  exact  only  if  c^  is  constant.  Fluid-property  laws  studied 
included:  first,  a  viscosity  proportional  to  temperature  and  an  arbitrary  but 
constant  Prandtl  number,  and  second,  the  "exact' *  property  laws  of  Cohen  (see 
Section  6.3).  The  assumption  of  viscosity  proportional  to  temperature  in  the 
boundary  layer  gives  a  value  of  C  =  p*i/peue  of  1.0,  since  density  is  in¬ 
versely  proportional  to  temperature  across  the  boundary  layer.  Recovery  factor 
versus  Mach  number  is  shown  in  figure  10a  for  the  various  laws.  Cohen's  re¬ 
lations  are  not  used  below  a  Mach  number  of  2.0  because  they  are  not  valid 
below  that  Mach  number  for  the  altitude  assumed,  50,000  ft.  They  are  not 
valid  for  W\ef  <  0.0152,  where  h^^  is  Cohen's  reference  enthalpy, 

2.119  x  10^  ft2/sec2.  The  assumptions  of  viscosity  proportional  to  temperature, 
C  =  1.0,  and  a  constant  Prandtl  number  of  0.72,  which  are  often  used  in  boundary, 
layer  work,  give  recovery  factors  that  are  within  3*5  percent  of  those  calcu¬ 
lated  with  the  more  nearly  exact  properties.  The  calculated  recovery  factors 
for  a  cone  are  identical  to  those  for  the  flat  plate  with  indentical  boundary 
conditions . 


Variation  of  the  calculated  recovery  factors  with  Prandtl  number  at  the 
wall  is  shown  in  figure  10b.  The  approximate  relation  that  is  often  assumed 
in  boundary-layer  theory 


rf  =  ~\f?r 


(7.11) 


is  also  shown. 
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7.5  Comparison  of  Present  Method  with  that  of  Flttgge-Lotz 

and  Blottner 

Recently  Flttgge-Lotz  and  Blottner  developed  an  implicit  finite-difference 
method  for  calculating  the  compressible  laminar  boundary  layer  (Reference  6). 
Fluid-property  relations  that  were  used  included  an  arbitrary  but  constant 
Prandtl  number  and  one  of  two  relations  for  viscosity.  Either  Sutherland's 
viscosity  law  or  a  viscosity  that  varied  linearly  with  temperature  was  used, 
and  examples  for  both  are  given  in  Reference  6.  Since,  as  of  now,  the  present 
method  is  programed  to  use  either  Cohen's  "exact"  relation  for  viscosity  or 
a  linear  variation  with  temperature,  results  of  the  two  methods  can  be  compared 
directly  only  for  the  examples  in  Reference  6  that  use  the  linear  viscosity 
law.  These  examples  consist  of  a  flat  plate  at  Mach  number  3.J  with  an  in¬ 
sulated  wall  and  with  a  heated  wall  such  that  T ,/T  ,  =  2.0.  Because  the 

vr  ad. 

boundary  conditions  are  independent  of  x,  the  flows  can  also  be  calculated 
by  similar-flow  methods.  Flttgge-Lotz  and  Blottner  compare  their  results  with 
the  similar  solution  of  Low  (Reference  l6).  The  results  of  the  present  method 
are  compared  with  the  same  method. 

The  shear  parameter  fj'  calculated  by  the  present  method  agrees  with 
that  calculated  by  Low's  method  within  the  accuracy  with  which  his  method  can 
be  used.  His  method  uses  parameters  presented  to  only  four  decimal  places. 

For  the  heated  plate  Low's  method  gives  fj'  «  0.35205  as  compared  to 
0.332057  obtained  by  the  present  method.  Flttgge-Lotz  and  Blottner  compare 
their  results  with  Low's  in  graphical  form,  and  it  appears  that  values  of 
fj'  agree  to  at  least  four  decimal  places  a  short  distance  downstream  of  the 
start  of  the  calculation.  At  the  start  the  agreement  between  the  two  methods 
is  not  so  close.  In  the  method  of  Flttgge-Lotz  and  Blottner  the  form  of  the 
equation  is  singular  at  x  =  0,  and  the  calculation  must  be  started  at  some 
point  downstream  of  x  =  0.  Furthermore,  the  initial  profile  must  be  specified 
at  this  starting  point.  Flttgge-Lotz  and  Blottner  used  the  similar  profile  of 
Low  to  get  started,  and  they  state  in  Reference  6  "Since  the  initial  profiles 
are  only  accurate  to  four  decimal  places,  there  are  initially  errors  in  the 
boundary -layer  calculations  and  characteristics.  These  errors  appear  to  de¬ 
crease  as  the  computations  proceed. ' ' 
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7.5  Comparison  of  Present  Method  with  that  of  Fltlgge-Lotz 

and  Blottner 

Recently  Fltlgge-Lotz  and  Blottner  developed  an  implicit  finite-difference 
method  for  calculating  the  compressible  laminar  boundary  layer  (Reference  6). 
Fluid-property  relations  that  were  used  included  an  arbitrary  but  constant 
Prandtl  number  and  one  of  two  relations  for  viscosity.  Either  Sutherland's 
viscosity  law  or  a  viscosity  that  varied  linearly  with  temperature  was  used, 
and  examples  for  both  are  given  in  Reference  6.  Since,  as  of  now,  the  present 
method  is  programed  to  use  either  Cohen's  "exact"  relation  for  viscosity  or 
a  linear  variation  with  temperature,  results  of  the  two  methods  can  be  compared 
directly  only  for  the  examples  in  Reference  6  that  use  the  linear  viscosity 
law.  These  examples  consist  of  a  flat  plate  at  Mach  number  3.0  with  an  in¬ 
sulated  wall  and  with  a  heated  wall  such  that  Tu/T^  =2.0.  Because  the 
boundary  conditions  are  independent  of  x,  the  flows  can  also  be  calculated 
by  similar-flow  methods.  Fltlgge-Lotz  and  Blottner  compare  their  results  with 
the  similar  solution  of  Low  (Reference  16).  The  results  of  the  present  method 
are  compared  with  the  same  method. 

The  shear  parameter  f^'  calculated  by  the  present  method  agrees  with 
that  calculated  by  Low's  method  within  the  accuracy  with  which  his  method  can 
be  used.  His  method  uses  parameters  presented  to  only  four  decimal  places. 

For  the  heated  plate  Low's  method  gives  f^'  =  0.33205  as  compared  to 
0.332057  obtained  by  the  present  method.  Fltlgge-Lotz  and  Blottner  compare 
their  results  with  Low's  in  graphical  form,  and  it  appears  that  values  of 
f^'  agree  to  at  least  four  decimal  places  a  short  distance  downstream  of  the 
start  of  the  calculation.  At  the  start  the  agreement  between  the  two  methods 
is  not  so  close.  In  the  method  of  Fltlgge-Lotz  and  Blottner  the  form  of  the 
equation  is  singular  at  x  =  0,  and  the  calculation  must  be  started  at  some 
point  downstream  of  x  =  0.  Furthermore,  the  initial  profile  must  be  specified 
at  this  starting  point.  Fltlgge-Lotz  and  Blottner  used  the  similar  profile  of 
Low  to  get  started,  and  they  state  in  Reference  6  "Since  the  initial  profiles 
are  only  accurate  to  four  decimal  places,  there  are  initially  errors  in  the 
boundary -layer  calculations  and  characteristics.  Hies e  errors  appear  to  de¬ 
crease  as  the  computations  proceed. ' ' 
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The  agreement  of  heat- transfer  values  calculated  by  the  present  method 
with  those  of  Low’s  similar  method  is  not  so  good  as  that  for  the  shear  para¬ 
meter.  For  the  heated  plate  at  Mach  3.0,  the  present  method  gave  a  value  of 
g^  of  -0. 29^16  as  compared  to  -0.29367  obtained  from  Low's  parameters. 
Since  Low  presents  his  data  to  four  decimal  places,  better  agreement  might  be 
expected.  The  disagreement  may  be  due  to  the  difference  in  the  recovery  fac¬ 
tor  calculated  and  vised  by  Low,  rf  =  0.8477,  and  that  calculated  by  the 
present  method,  0.84746.  Flt!gge-Lotz  and  Blottner  also  compare  their  heat- 
transfer  values  with  Low's  values  calculated  for  the  same  flow.  Close  obser¬ 
vation  of  their  graphs  at  a  point  downstream  of  the  initial  x-station  appears 
to  show  that  their  results  are  around  0.2  percent  higher  than  Low's  (see 
figures  6(b)  and  7(b)  of  Reference  6).  This  is*  about  the  same  percentage 
difference  that  appears  between  the  present  method  and  Low's,  that  is. 

Present  method  g.' 

- —  =  1.002 

Low's  g^ 

Velocity  and  enthalpy  profiles  as  calculated  by  the  three  methods  —  Low, 
Flttgge-Lotz  and  Blottner,  and  present  —  agree  to  three  decimal  places  at  the 
one  station,  x  =  1.10,  where  data  exist  so  that  such  a  comparison  could  be 
made.  This  is  within  the  accuracy  with  which  the  profiles  are  presented  in 
Reference  6. 


7.6  Comparison  of  Heat  Transfer  Calculated  by  the  Present 
Method  with  that  Measured  on  a  Circular  Cylinder 

Few  experimental  data  on  compressible  boundary  layers  sure  available  for 
comparison  with  the  prediction  of  theories.  Exceptions  sure  a  few  heat-trans¬ 
fer  data.  Reference  17  presents  experimental  data  on  heat  transfer  measured 
on  circular  cylinders  with  their  axes  normal  to  the  flow  at  several  Mach 
numbers.  Comparison  of  the  measured  heat  transfer  with  that  calculated  by  the 
present  method  is  shown  in  figure  11.  Rie  external  velocity  and  onset  flow 
conditions  used  in  the  calculation  were  obtained  from  Reference  17.  The  ex¬ 
perimental  temperatures  were  so  low  that  the  relations  for  fluid  properties 
given  by  Cohen  (see  Section  6.3)  are  invalid.  Therefore,  the  simple  assunp- 
tions  of  constant  Prandtl  number  (0.72)  and  a  viscosity  proportional  to 
temperature  through  the  boundary  layer  were  assumed. 
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As  discussed  by  the  authors  of  Reference  17,  the  magnitude  of  the  experi¬ 
mental  error  is  not  known;  but  it  may  be  as  great  as  the  difference  between 
the  calculated  and  measured  results  in  figure  11.  They  point  out  that  the 
experimental  data  at  M  =  4.15  are  known  to  include  experimental  error  not 
present  at  the  other  Mach  numbers,  and,  furthermore,  that  the  error  in 
instrumentation  tends  to  increase  the  indicated  heat  transfer  rearward  of  *+0®. 
They  also  state  that  "for  a  ^  90°  the  accuracy  of  all  the  data  becomes 
questionable  because  of  the  very  small  convective  heat  rates  compared  with 
radiation  and  other  sources  of  error."  To  avoid  confusion  with  experimental 
points,  calculated  points  are  not  indicated  in  the  figure,  but  the  step  size 
was  5°  except  near  the  stagnation  point, where  step  lengths  of  1®  were  used, 
analler  steps  were  used  here  in  order  to  converge  on  the  heat-transfer  rate 
at  the  stagnation  point.  Values  calculated  are  believed  to  be  accurate  to 
three  decimal  places.  Experimental  data  were  obtained  at  gw  —  1.0,  0.8, 
and  0.7«  The  calculated  curves  shown  are  for  =  0.8,  but  calculations  at 
the  other  two  enthalpies  were  very  similar.  A  gw  of  0.70  gives  values  of 
fs«/ (Sg^  ~  gw)]/Cgw/(gad  -  that  are  from  0.7  to  1.2  percent  higher 

than  those  shown. 

7-7  Boundary  Layer  in  an  Axis ymme trie  Convergent-Divergent 

Rocket  Nozzle 

As  an  example  of  internal  axisymmetric  flow  the  boundary  layer  inside 
a  rocket  nozzle  was  calculated.  The  flow  involves  both  large  pressure 
gradients  and  heat  transfer.  It  is  similar  to  the  flow  calculated  by  C.  B. 
Cohen  and  E.  Reshotko  in  Reference  18,  but  an  exact  comparison  between  the  two 
calculations  can  not  be  made  because  there  is  insufficient  information  on  the 
flow  properties  and  the  nozzle  dimensions  of  Reference  18.  Radius,  surface 
distance,  and  the  velocity  and  Mach  number  at  the  center  of  the  nozzle  are 
plotted  in  figure  12  against  the  axial  distance.  These  local  flow  conditions 
were  obtained  by  means  of  simple  one-dimensional  flow  relations.  Stagnation 
pressure  is  assumed  to  be  500  pounds  per  square  inch  and  stagnation  tempera¬ 
ture  as  4000®R,  which  leads  to  a  total  enthalpy  Hg  =  27.246  x  10  ft^/sec^. 
The  nozzle  wall  is  assumed  cooled  to  a  uniform  temperature  of  800*R,  which 
corresponds  to  an  enthalpy  ratio  of  gw  =  0. 17625. 
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Results  of  the  calculation  are  presented  in  figure  13.  The  skin-friction 
coefficient  c^  ,  defined  by  equation  (6.115)#  and  the  ratio  of  Stanton  number 
to  skin  friction  coefficient#  defined  by  equations  (6.120)  and  (6.121)#  are 
shown.  Calculations  denoted  by  1 'exact  properties  1 '  were  obtained  by  using 
the  relations  of  N.B. Cohen  described  in  Section  6.3*  A  calculation  was  also 
made  by  using  the  simple  assumptions  of  constant  Prandtl  number  and  viscosity 
proportional  to  local  enthalpy.  Because  of  the  high  enthalpies  involved,  a 
Prandtl  number  of  0.78  was  chosen.  Figure  13  shows  that  these  simple  assump¬ 
tions  for  fluid  properties  lead  to  results  that  seriously  overestimate  both 
skin  friction  and  heat  transfer  at  the  higher  Mach  number.  At  the  lower  Mach 
numbers,  the  effect  is  reversed#  that  is#  the  simple  assumptions  gave  lower 
skin  friction  and  heat  transfer  than  those  calculated  for  exact  properties. 
This  study  and  those  in  previous  sections  indicate  that  the  exact-property 
relations  should  be  used  if  the  local  enthalpy  exceeds  4  x  10  ft2/sec2. 

Stagnation-point  flow  is  assumed  at  the  start  of  the  flow  (x  =  0). 

Because  ufi  =  0,  here,  the  skin  friction  coefficient  cf  is  singular  and, 
of  course,  can  not  be  presented  at  the  stagnation  point  in  figure  13-  The 
shear  parameter  fy  and  also  the  gradient  of  g  at  the  wall,  gj  ,  are 
finite  here  and  are 

fj'=  0.38084  j  =  0.16523  for  exact  properties 

f”  =  0.73507  J  =  0.38177  for  C  =  1.0,  Pr  =  O.78 


7.8  Boundary  Layer  on  a  Spherically  Blunted  Cone  at  Mach 

Numbers  of  3.0  and  9*0 

Boundary  layers  were  calculated  on  a  reentry- type  body  that  consists  of 
a  15 0  half-angle  cone  blunted  with  a  spherical  nose.  The  conical  and  spherical 
surfaces  join  in  such  a  way  that  there  is  no  discontinuity  in  body  slope.  The 
velocity  distribution  outside  the  boundary  layer  was  calculated  by  a  computer 
program  that  consists  of  (l)  a  solution  of  the  blunt-body  problem  from  the 
stagnation  point  to  a  point  past  the  sonic  point  and  (2)  a  solution  by  the 
method  of  characteristics  in  the  supersonic  flow  region.  The  method  of 
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characteristics  Is  started  Just  aft  of  the  sonic  point  by  using  local  veloci¬ 
ties  and  pressures  given  by  the  blunt-body  solution  as  an  initial  boundary 
condition.  Variation  of  the  inviscid  velocity  and  Mach  number  with  axial 
distance  is  shown  in  figure  14.  Surface  distance  and  local  radius  are  also 
shown.  Though  the  velocity  varies  smoothly  in  going  from  the  blunt-body 
solution  to  the  characteristics  solution,  there  is  a  large  change  in  velocity 
gradient,  dufi/dx,  just  aft  of  the  point  where  the  change  in  methods  occurs, 
x  =  0.9.  This  results  in  large  changes  in  the  parameter  P  =  (x/ue)(due/dx) 
just  aft  of  this  point,  which  produces  large  changes  in  the  calculated  values, 
as  will  be  seen  below. 


Boundary  layers  on  the  body  were  calculated  at  the  two  speeds,  M  =  3.0 
amd  9*0,  at  three  altitudes,  50,000,  100,000  and  250,000  ft.,  for  three  body 
temperatures.  Results  for  M  =  3.0  are  shown  in  figure  15  (a,  b,  c).  Results 
are  independent  of  altitude  to  three  decimal  places;  therefore, only  the  50,000- 
ft  values  are  shown.  The  reason  that  the  effect  of  altitude  is  found  to  be 
negligible  is  that  there  is  only  a  small  change  in  total  enthalpy  with  altitude. 


Total  enthalpy  is  given  in  Table  VII  for  the  three  altitudes.  Also,  the 
reference  velocity  a^.  that  is  used  to  nondimensionalize  the  velocity  ug  is 
given.  The  velocity  a*  is  the  velocity  of  sound  at  free-stream  stagnation 
conditions.  The  total  enthalpy  of  about  6.6  x  10  ft  /sec  at  the  lower 
altitudes  corresponds  to  a  total  temperature  of  about  1090 °R.  The  three  body 


temperatures  presented  are:  -  0.357,  which  gives  a  body  temperature  equal 

to  about  the  static  temperature  of  the  onset  flow;  gw  =  0.4-58,  which  gives  a 
wall  temperature  equal  to  about  ^0°F)  and  an  adiabatic  wall,  which  is  calculated 
by  setting  =  0.  Variation  of  the  skin-friction  coefficient  cf  and  ratio 
of  Stanton  number  to  the  skin  friction  coefficient  St/c±,  with  surface  distance 
are  shown  for  the  flows  when  is  specified.  When  an  adiabatic  wall  (g^  =  0) 

is  specified,  c^  and  are  shown.  To  further  investigate  the  effect  of 
fluid-property  laws  on  the  solution,  the  flows  were  calculated  twice,  once 
with  the  "exact  fluid  properties"  of  Cohen  (Section  6.5),  and  again  with  the 
simple  assumptions  of  constant  Prandtl  nunber  and  a  linear  viscosity  law.  It 


is  seen  that  the  simple  assumptions  greatly  overestimate  the  skin  friction  and 
heat  transfer  (figure  15(a,  b))  and,  in  the  adiabatic  wall  case,  the  surface 
temperature  (figure  15c).  For  the  results  with  exact  properties  the  figures 


show  a  "bump"  in  c^  just  aft  of  x  =  0.90.  This  is  at  the  point,  mentioned 
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in  the  last  paragraph,  where  the  method  of  calculating  the  inviscid  flow  re¬ 
sults  in  a  large  gradient  in  ug.  This  gradient  causes  a  large  variation  in 
the  local  enthalpy  h  in  this  region,  which  produces  a  large  variation  in 
Prandtl  number  (see  figure  2).  Even  though  there  is  a  large  change  in  local 
flow  conditions  it  is  believed  that  the  calculated  values  are  accurate  to  three 
decimal  places,  because  in  accuracy  studies  where  the  present  method  was  vised 
to  calculate  flows  with  discontinuities  in  P,  and  also  in  gw  and  fy,  such 
accuracy  was  obtained  downstream  of  the  discontinuity. 

Results  for  M  =  9*0,  shown  in  figure  l6(a,  b,  c),  have  a  variation  similar 
to  those  at  the  lower  speed.  As  was  also  the  case  for  the  flow  at  Mach  number 
3.0,  altitude  has  a  small  effect  on  the  calculated  values.  A  change  from 
50,000  ft  to  250,000  ft  causes  at  most  a  change  of  4  in  the  third  significant 
figure  of  the  values  calculated.  Total  enthalpy  and  the  reference  velocity  vised 
in  the  calculations  are  given  in  Table  VII.  The  total  enthalpy  of  about 
4.0  x  10  ft  /sec  corresponds  to  a  total  temperature  of  about  6700°R.  Since 
some  dissociation  would  occur  at  this  high  a  temper ature,  there  would  be  some 
change  in  the  fluid  properties.  But  the  magnitude  of  this  change  and  its  effect 
on  the  solution  are  left  for  a  later  report.  It  will  be  simple  to  study  the 
effect  of  equilibrium  dissociation  bymeans  of  the  present  method  of  solution, 
since  all  that  will  be  necessary  is  to  replace  the  functions  for  Prandtl  number 
given  in  Section  6.3  by  similar  functions  for  an  effective  Prandtl  number  that 
includes  dissociation.  Such  a  Prandtl  number  is  presented  in  Reference  J.  The 
three  body  temperatures  in  figure  16  are:  g^  =  0.083124,  which  gives  a  wall 
temperature  of  about  560°^  or  100“F;  gw  =  0.40,  which  gives  a  wall  temperature 
of  around  2000 °R;  and  the  adiabatic  wall,  which  is  given  by  =  0.  Again 
the  simpler  fluid-property  laws  —  those  of  a  viscosity  that  varies  linearly 
with  temperature  and  a  constant  Prandtl  number  —  give  results  that  differ 
greatly  from  those  obtained  with  the  "exact  properties".  Values  of  St/cf 
are  indeterminate  at  the  stagnation  point,  but  values  of  f^'  ,  gtf,  and  g^ 
are  finite  and  are  given  in  Table  VIII. 
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TABLE  VII 


ALTITUDE,  TOTAL  ENTHALPY,  AND  REFERENCE  VELOCITY  USED  IN 
CALCULATION  OF  FLOWS  OVER  THE  SPHERICALLY  BLUNTED  CONE 


Altitude 

Static 

Mach  3 

Mach 

9 

ft 

Temperature 

H  x  10-T 
e 

a ,  =  u 
t  » 

H  x  10"7 
e 

at  =  u» 

®R 

2 /  2 

ft  /sec 

ft /sec 

. ,  2  /  2 

ft  /sec 

ft/sec 

50,000 

390.0 

0.65523 

1620.5 

4.03180 

4016.5 

100,000 

392.4 

O.65893 

I625.O 

4.05366 

4027.4 

250,000 

451.8 

0.75894 

1743.6 

4.66769 

4321.5 

TABLE  VIII 


STAGNATION  POI17T  VALUES  OF  fj1,  g^,  AND 


FOR  THE 


SPHERICALLY  BLUNTED  CONE 


Moo 

Altitude 

ft 

fll 

w 

Exact 

Properties 

Pr  =  0.72 
C  =  1.0 

Exact 

Properties 

Pr  =  0.72 
C  =  1.0 

Exact 

Properties 

Pr  =  0.72 
C  =  1.0 

3.0 

50,000 

1.3119 

1.0 

1.0 

* 

0.0 

* 

0.0 

3.0 

50,000 

O.985I 

0.9828 

0.3570* 

0.3570* 

0.4205 

0.4076 

3.0 

50,000 

1.1365 

1.0369 

0.4580* 

O.4558 

O.3878 

0.3475 

9.0 

50,000 

0.9261 

1.3119 

1.0 

1.0 

* 

0.0 

* 

0.0 

9.0 

50,000 

0.3155 

0.8275 

0.08312* 

0.08312* 

0.1883 

0.5603 

9.0 

50,000 

0.6091 

0.400* 

0.400* 

0.2117 

0.3822 

9.0 

250,000 

0.9261 

1.3119 

1.0 

1.0 

0.0* 

0.0* 

9.0 

250,000 

0.3185 

0.8275 

0.08312* 

0.08312* 

0.1854 

0.5603 

9.0 

250,000 

0.6113 

l.or6o 

0.400* 

0.400* 

0.2162 

0.3822 

Input  values 
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7.9  Transverse  Curvature  Effect  in  Axially  Symmetric  Flow 


In  all  the  examples  presented  so  far,  the  transverse  curvature  terms  have 

been  neglected.  These  terms  become  of  importance  when  the  boundary  layer  has 

a  thickness  of  about  the  same  magnitude  as  the  body  radius.  There  are  two  of 

these  terms:  the  first  term.,  which  includes  T,  on  the  right-hand  side  of 

both  the  momentum  equation  (6.35)  and  the  energy  equation  (6.46)  and  the  term 
p  +  1 

that  contains  N  =  — j —  +  R  in  both  equations.  The  transverse  curvature 
parameter  T  is  defined  by  (6.32)  and  the  radius  parameter  R  by  (6.35b). 


Transverse  curvature  effects  were  investigated  by  Probstein  and  Elliott, 
Reference  19.  The  study  was  restricted  to  similar  flow  and  essentially  ex¬ 
panded  the  stream  function  f  and  the  enthalpy  function  g  in  a  power 
series  in  |,  where 


* 

6  cos  a 


(7.12) 


That  is, 


f(t»  n)  =  fj(n)  +  1  ^(n)  +  t2  f2(n)  ••• 


(7.15) 


and  there  is  a  similar  function  for  enthalpy.  The  expansion  parameter  was 
assumed  to  be  small  with  respect  to  unity.  Substituting  their  asymptotic 
expansions  into  the  momentum  and  energy  equations,  and  equating  to  zero  all 
terms  with  the  same  power  in  5,  Probstein  and  ELliot  obtained  a  double 
infinity  of  ordinary  differential  equations .  They  found  all  these  equations 
except  the  zeroth  order  momentum  equation  to  be  linear,  because  a  Pr  =  1.0 
and  a  viscosity  proportional  to  temperature  had  been  assumed.  They  obtained 
solutions  of  the  zeroth-  and  first-order  equations.  They  forced  the  equations 
to  be  similar,  that  is,  independent  of  x,  by  assuming  i  to  be  constant. 

As  discussed  by  them,  this  assumption  requires  that  the  boundary-layer  thick¬ 
ness  be  directly  proportional  to  the  body  radius.  They  present  solutions  for 
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only  two  types  of  bodies  —  a  cone  and  a  cylinder  with  axes  aligned  with  the 

* 

flow.  Their  assumption  of  similarity  —  that  8  varies  as  the  body  radius  — 
is  a  poor  approximation  for  these  two  bodies,  for  in  reality  the  thickness 
grows  nearly  parabolically  with  x,  as  V3T. 


The  present  method  was  used  to  calculate  the  boundary  layer  on  the  same 
two  types  of  bodies.  Assumptions  on  the  fluid  properties  are  the  same  as 
those  used  by  Probstein  and  Elliott.  Results  are  presented  in  terms  of  the 
parameters  used  by  them, 

- 1 jJSjlO  (7.1*0 

(fw'^T  =  0 

for  the  cylinder,  and 

f  1 1  _  (f  1 1 ) 

taa  a-\/lie  f.X-O  (7.15) 

v  w  'T  =  0 


for  the  cone,  where  Re  =  ugx/ve  >  a  i8  the  angle  between  the  normal  to  sur¬ 
face  and  the  radius,  but  is  also  the  half  angle  of  the  cone  here;  and  the  sub¬ 
script  T  =  0  denotes  the  solution  with  no  transverse  curvature  effect.  For 
the  cylinder  the  radius  parameter  R  is  zero,  and  for  the  cone  it  is  given 
exactly  by  (6.35b) 


R 


x  drp 

rn  3*T 


1  + 


*/ 


dTJ 


(6.33b) 


since  cos  a  is  constant.  Calculations  for  a  large  range  in  Reynolds  numbers, 

•y 

up  to  10  ,  at  Mach  numbers  from  1.0  to  10.0,  wall  enthalpy  ratios  hw/he  from 
1.0  to  10.0,  and,  for  the  case  of  cone  flow,  a  of  from  5*  to  20°  showed  that 
the  parameters  (7*1*0  and  (7*15)  asymptotically  approached  a  constant  as  the 
solution  proceeded  downstream.  It  was  found  that  the  asymptotes  depended 
mostly  only  on  Mach  number  and  body  temperature  and  were  nearly  independent 
of  Reynolds  number,  body  radius,  and,  in  the  case  of  the  cone,  the  cone  angle. 
Sample  solutions  are  shown  in  figures  17  and  18  for  a  cylinder  and  cone, 
respectively.  Probstein  and  Elliott's  solutions  for  the  same  flows  are  indi¬ 
cated  in  the  figures  and  give  values  of  the  parameters  plotted  that  are 
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Figure  17.-  Effect  of  transverse  curvature  on  shear  parameter  for  flow  over  a  circular  cylinder 
with  axis  aligned  with  the  flow.  Me=  5.0;  hw/he=  1.0;  Pr  =  1.0. 


Figure  18  -  Effect  of  transverse  curvature  on  shear  parameter  for  flow  over  a  cone  with  axis 
aligned  with  the  flow.  Me  =  5.0;  hw/he=  5.0;  Pr  =  1.0;  a  =  10°;  Re  =  104/ft. 
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independent  of  x.  Because  of  the  similarity  requirement  assumed  by  PrOb- 
stein  and  Elliott,  the  disagreement  between  the  two  methods  is  not  sur¬ 
prising.  In  the  calculations  by  the  present  method,  solutions  were  required 
to  converge  to  four  significant  digits. 


Effects  of  Mach  number  and  wall  temperature  on  the  asymptotic  values  of 
the  shear  parameter  are  shown  in  figure  19  for  the  cone.  The  Probstein- 
Elliott  values  are  also  shown.  Although  the  values  given  by  the  two  methods 
of  solution  differ,  the  effects  of  Mach  number  and  wall  temperature  are 
similar.  Since  a  Prandtl  number  of  1.0  is  assumed  for  the  flows  studied, 
Reynolds  analogy  holds  exactly,  and  therefore  the  effect  of  transverse  curva¬ 
ture  on  heat  transfer  is  identical  to  its  effect  on  wall  shear,  that  is. 


f » »  _  (r  t » ) 

w  uw  ;T  =  0 

(f » t ) 

v  w  'T  =  0 


(7.16) 


As  was  pointed  out  by  Probstein  and  Elliott,  transverse  curvature  behaves  like 
an  external  favorable  pressure  gradient,  that  is,  it  increases  both  local  shear 
and  heat  transfer  and  tends  to  delay  both  separation  and  transition. 


7.10  Effect  of  Discontinuities  in  Surface  Temperature  and 
in  Wall  Mass  Transfer  on  Wall  Shear  and  Heat  Transfer 


In  Reference  2  the  authors  investigated  solutions  of  the  incompressible 
boundary  layer  in  the  region  where  a  discontinuity  in  the  pressure  gradient 
parameter  P  occurred.  The  purpose  of  the  investigation  was  to  learn  how 
rapidly  the  boundary  layer  adjusts  to  local  changes  in  pressure  outside  the 
boundary  layer.  The  answer  to  this  question  is  of  great  importance  in  appli¬ 
cation  of  local -similarity  methods  that  give  approximate  solutions  of  the 
boundary-layer  equations.  In  local-similarity  methods  solutions  are  obtained 
by  a  step-by-step  procedure,  in  which  each  x-wise  segment  of  the  flow  is 
approximated  by  one  of  the  family  of  similar  flows.  The  method  assumes  that 
the  local  boundary-layer  profiles  are  functions  of  local  boundary  conditions 
only.  For  the  method  to  be  accurate  the  boundary  layer  must  adjust  to  the 
local  boundary  conditions  very  rapidly. 
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In  the  Investigation  of  the  incompressible  boundary  layer  in  the  region 
of  a  discontinuity  in  P,  it  was  found  that  this  adjustment  to  local  boundary 
conditions  was  surprisingly  slow  (Reference  2).  Study  of  the  same  adjust¬ 
ment  for  compressible  flow  by  means  of  the  present  method  shows  that  it  is 
nearly  identical  to  that  found  earlier  for  incompressible  flow.  Also,  the 
present  method  was  used  to  find  how  rapidly  the  boundary  layer  adjusted  to 
discontinuities  both  in  wall  temperature  and  in  mass  transfer  at  the  surface 
for  supersonic  flow.  These  effects  are  presented  in  this  Section. 

First,  the  effect  of  discontinuities  in  wall  temperature  gtf  on  a  flat 
plate  at  Mach  number  3.0  are  presented.  Prandtl  number  is  assumed  to  be  0.72. 
The  wall  temperature  is  constant  forward  of  x  *  1.0  and  then  abruptly 
changes  to  another  value  downstream  of  x  -  1.0.  Three  values  of  the  wall 
temperature  were  studied:  hy  =  h^,  giving  a  value  of  gtf  *  0.90134; 

hw  -  2  h^,  giving  gy  =  1.80268;  and  hy  =  hg,  giving  gw  =  0.35714.  In 

figure  20a  calculated  values  of  g^  are  shown  for  the  flows  when  gy  is 
initially  the  adiabatic  value  and  then  abruptly  changes  to  the  higher  or 
lower  temperature  downstream  of  x  »  1.0.  For  comparison,  the  calculated 
values  are  also  shown  for  the  case  when  the  plate  is  adiabatic  for  the  whole 
length  of  flow.  The  latter  is  a  similar  flow  and  therefore  independent  of  x. 

A  positive  value  of  g^  means  that  heat  1b  being  transferred  to  the  surface; 
a  negative  value  means  that  heat  is  being  transferred  from  the  surface  to  the 

air.  The  figure  shows  that  abruptly  changes  just  aft  of  the  discorf- 

tinuity  in  gw  and  slowly  approaches  an  asymptotic  value  as  the  flow  proceeds 
downstream.  This  asymptotic  value,  which  is  indicated  in  the  figure,  is  the 
same  as  that  which  would  be  obtained  by  a  similar  solution  using  the  local 
boundary  conditions  downstream  of  the  discontinuity.  It  is  seen  that  the  cal¬ 
culated  value  of  gj  overshoots  the  asymptotic  value  just  downstream  of  the 
discontinuity  and  slowly  approaches  the  value  from  the  opposite  side  from  which 
it  started.  That  is,  in  figure  20a  for  the  case  whbn  gw  =  0.35714  downstream, 
the  calculated  value  of  gj  at  first  greatly  exceeds  its  asymptotic  value  and 
then  approaches  it  from  above  as  the  calculation  proceeds  downstream.  Step 
size  in  x  used  in  the  calculation  was  selected  in  such  a  way  as  to  give  cal¬ 
culated  values  accurate  to  three  decimal  places  downstream  of  the  discontinuity. 
Local  similarity  methods  would  give  a  constant  value  of  for  the  whole 
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length  of  flow  downstream  of  the  discontinuity ,  and,  furthermore,  this  value 
would  be  the  asymptotic  value  indicated  in  the  figure. 

Similar  plots  for  flows  where  gy  is  initially  0.35714  and  1.80268 
and  then  abruptly  changes  to  the  other  two  values  are  shown  in  figures  20b 
and  20c,  respectively.  The  adjustment  of  g^  to  the  local  boundary  condi¬ 
tions  is  like  that  shown  in  figure  20a.  As  was  the  case  in  the  adjustment 
of  fj'  to  local  changes  in  P,  the  adjustment  of  g^  to  the  local  value  of 
gw  is  surprisingly  slow.  Since  flat-plate  flow,  P  =  0,  is  assumed  for  the 
flows  shown  in  figure  20  (a,  b,  and  c),  the  term  containing  pg/p  in  the 
momentum  equation  (6. 37)  drops  out,  and  the  calculated  values  of  f,  f •,  and 
f  *  *  are  independent  of  gy. 

Brthalpy  profiles  for  two  of  the  flows  with  discontinuity  in  gw  are 
shown  in  figure  21.  It  is  seen  that  the  enthalpy  values  adjust  to  the  local 
boundary  conditions  more  rapidly  near  the  wall  than  they  do  at  the  outer  edge 
of  the  boundary  layer.  Asymptotic  profiles  for  the  two  flows  are  also  shown. 

Discontinuities  in  mass  transfer  at  the  surface  were  also  investigated. 
Two  external  flows  were  studied.  First  flat-plate  flow  similar  to  that  des¬ 
cribed  in  the  previous  paragraphs . is  assumed:  P  =  0.0,  Mg  =  3*0,  Pr  =  0.72. 

It  is  assumed  that  air  is  blowing  out  of  the  surf  we  forward  of  x  =  1.0  in 
such  a  way  that  f  is  equal  to  a  constant.  Aft  of  x  =  1.0  the  surface 
is  impermeable.  Although  there  is  a  discontinuity  in  the  suction  velocity  vw 
the  stream  function  f  is  a  continuous  function  since  it  is  an  integral  of 
vw  over  x  (see  equation  6.1 6).  Two  different  blowing  quantities  are  assumed 
forward  of  x  =  1.0:  f  »  -0.5  and  -0.1.  The  resulting  variation  of  f 
downstream  of  the  discontinuity  in  vy  is  shown  in  figure  22a.  Effects  of 
the  mass  transfer  and  its  discontinuity  on  the  calculated  values  of  fj'  and 
gy  are  shown  in  figures  22b  and  22c,  respectively,  for  the  case  of  an  in¬ 
sulated  plate.  It  is  seen  that  blowing  decreases  both  i£*  and  g^. 

The  effect  of  blowing  cold  air  at  the  front  of  a  nearly  insulated  plate 
is  shown  in  figure  23.  The  same  blowing  quantities  used  in  the  previous  para¬ 
graph  arc  assumed  (figure  22a).  The  air  blown  from  the  surface  is  assumed  to 
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have  the  sane  temperature  as  the  free  stream,  h  =  h  ,  giving  a  value  of 
gy  =  0.3571^  forward  of  x  =  1.0.  Aft  of  x  =  1.0,  the  surface  is  solid 
and  has  a  temperature  close  to  the  adiabatic  temperature.  Effects  of  the 
blowing  and  its  discontinuity  on  are  shown  in  figure  23b.  For  com¬ 
parison  values  of  for  the  case  of  no  blowing  (fy  =  0)  are  also  plotted. 
Since  flat-plate  flow  (P  =  0)  is  assumed,  the  term  containing  p^p  in  the 
momentum  equation  disappears,  and  the  calculated  values  of  f,  f',  and  f'1 
are  independent  of  gy.  Therefore  the  variation  of  fj'  with  x  for  this 
case  of  blowing  cold  air  on  a  hot  plate  is  the  same  as  that  shown  in  figure 
22b  for  the  case  of  blowing  on  an  insulated  plate. 

The  opposite  effect,  that  of  blowing  hot  air  at  the  front  of  a  cold  plate, 
is  shown  in  figure  2k.  Three  different  blowing  quantities  are  assumed  forward 
of  x  =  1.0:  f ^  =  —  0.1,  —  0.5,  and  —  1.0.  The  first  two  of  these  are 
identical  to  those  used  on  the  plates  described  in  the  two  previous  paragraphs 
and  shown  in  figure  22a.  The  blowing  air  is  assumed  to  have  a  temperature 
nearly  equal  to  the  adiabatic  temperature,  gw  =  1.0.  Downstream  of  blowing 
(x  >  1.0)  the  surface  has  a  temperature  equal  to  the  static  temperature  out¬ 
side  the  boundary  layer.  Again  the  variation  of  gj  with  distance  is  shown. 
The  effect  of  blowing  on  the  velocity  and  enthalpy  profiles  for  this  flow  is 
shown  in  figure  25.  The  adjustment  of  the  profiles  to  the  discontinuities  in 
mass  transfer  and  temperature  is  typical  of  that  observed  for  all  flows  where 
such  discontinuities  occurred. 

The  effect  of  blowing  air  out  of  the  surface  at  the  front  of  a  body  is 
applicable  to  some  forms  of  ablation  cooling.  In  practice,  when  ablation 
cooling  is  used,  the  body  is  usually  blunt  and  thus  the  flow  is  initially  like 
that  of  a  stagnation-point  flow  rather  than  the  flat-plate  flow  that  was 
assumed  in  the  preceding  paragraphs.  Therefore  the  effect  of  blowing  has  also 
been  studied  on  a  body  that  is  initially  blunt.  The  flow  is  assumed  to  be 
stagnation-point  flow  (P  =  1.0)  for  x  i  1.0  and  then  to  change  abruptly  to 
a  constant  pressure  flow  (P  =  0.0)  for  x  >  1.0.  This  assumption  of  an 
abrupt  change  from  stagnation-point  flow  to  constant-pressure  flow  would  not 
be  a  bad  approximation  to  the  flow  over  the  reentry  type  body,  a  spherically 
blunted  cone,  studied  in  Section  7*8  and  shown  in  figure  lk.  The  mass-trans¬ 
fer  distributions  are  assumed  to  be  the  same  as  those  studied  in  the  flat- 
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plate  flow.  Also  the  onset  flow  is  the  same.  Effect  of  the  discontinuities 
in  P  and  vy  on  an  insulated  body  is  shown  in  figure  26.  Variation  of 
the  shear  parameter  f^'  with  distance  is  shown  in  figure  26c  and  of  the 
mil  enthalpy  gw  with  distance  in  figure  26d.  Rote  that  even  though  the 
pressure  gradient  parameter  P,  the  suction  parameter  fy,  and  the  heat- 
transfer  parameter  g£  are  constants  in  the  stagnation-point  flow  region, 
x  <  1.0,  the  wall  enthalpy,  and  thus  f^',  change  because  of  the  variation  in 
u .  The  effect  of  blowing  cold  air  initially  out  of  the  blunt  body  is  shown 

C 

in  figure  27.  Aft  of  the  discontinuities  in  both  F  and  vy,  the  surface  is 
assumed  to  have  a  temperature  equal  to  the  total  temperature  of  the  flow. 
Adjustment  of  the  velocity  and  enthalpy  profile  to  the  change  in  P,  vy,  and 
gy  is  shown  in  figure  28  for  one  of  the  mass -transfer  distributions  of 
figure  27. 

In  conclusion,  the  study  of  discontinuities  in  the  boundary  conditions  — 
P,  g^ ,  and  fw  —  shows  surprisingly  slow  adjustment  of  the  boundary  layer  to 
these  local  conditions.  It  appears  that  local  similarity  methods  give 
erroneous  results  for  flows  where  rapid  variation  in  the  external  or  wall 
boundary  condition  occurs. 

7.11  Concluding  Statement 

The  large  nuriber  of  calculations  and  their  comparison  with  other  methods 
and  experiment  establish  the  fact  that  the  method  is  rapid,  highly  accurate, 
and  powerful.  About  0.25  minutes  of  computing  time  on  the  IBM  7090  computer 
are  required  to  calculate  the  velocity  and  enthalpy  profiles  at  one  station 
accurate  to  four  significant  digit.  Since  a  typical  flow  requires  from  20 
to  30  stations,  computer  time  averages  between  5  aoi  7  minutes  for  solving  a 
given  problem.  The  method  appears  capable  of  solving  any  flow  problem  for 
which  the  boundary-layer  equations  themselves  remain  valid.  It  is  now  being 
extended  to  handle  flows  where  dissociation  (both  equilibriun  and  nonequi¬ 
librium  )  occurs.  Also  effects  of  both  catalytic  and  noncatalytlc  walls  will 
be  Included. 
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APPENDIX  A 


DEVELOPMENT  OP  THE  BOUNDARY -LAY®  EQUATIONS  FOR  COMPRESSIBLE 
LAMINAR  FLOW  INCLUDING  THE  EFFECT  OF  TRANSVERSE  CURVATURE 

In  this  Appendix  the  boundary-layer  equations  are  developed  for  axially 
symmetric  flow  for  the  coordinate  system  of  Section  6.1.  The  effects  of 
transverse  curvature  are  included. 

Continuity  Equation 

The  continuity  equation  is 

||  +  div  (p  V)  =  0  (Al) 

The  first  term  is  of  course  zero  for  steady  flow. 

Before  writing  (Al)  in  the  x,y  coordinate  system,  consider  the  general 
orthogonal  coordinates .  Let 

X,  Y,  Z,  09  any  Cartesian  coordinates 
aV  a2*  (xya°  any  ^"thosonal  coordinates 

Let  the  elements  of  length  in  the  direction  of  increasing  ct^,  a^,  and 
be  hjCt^,  hgOtg,  and  h^a^,  respectively,  so  that 

(ds)2  =  hL2(da1)2  +  h22(dfl2)2  +  h?2(d a^)2  (A2) 

The  quantities  h^,  hg,  and  h^  are  defined  by 
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Now  consider  the  coordinate  system  of  figure  1.  It  is  redrawn  in  figure  Al, 
to  show  that 


dotj.  -  dx 

dflfg  ■  dy 
dcij  =  dtp 

where  here  <p  is  the  azimuthal  angle.  If  the  second  order  effect  of  longi¬ 
tudinal  curvature  were  to  be  retained  as  well  as  that  of  transverse  curvature, 
=  dx/(l  +  ky)  where  k  is  the  local  curvature,  that  is,  the  reciprocal 
of  the  local  radius  of  curvature.  For  most  flows  of  interest  in  aerodynamics 
the  effect  of  longitudinal  curvature  is  usually  much  smaller  that  that  of 
transverse  curvature,  and  it  will  be  neglected  in  the  present  study.  Since 
the  problem  is  axlsymmetric,  all  quantities  are  independent  of  q>.  Now  from 
(A3) 

■  hg  =  1 
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Prom  vector  analysis  and  the  relation 

V  =Tu  +  J  v  (A5) 

the  continuity  equation  can  be  written 

§*+div  (PV)a|  +  i^(rpu)  +  ^(rpv)}=0  (A6) 

*;  o 

Momentum  Equation 

The  Navier- Stokes  equation  of  motion  in  vector  notation  is  (from  page  50 
of  Howarth,  Reference  5)i 

p  —  V  X  curl  V  +  i  grad  V  2  j«  p  F  —  grad  P  +  j  gr®d  (uA)  +  grad  (V  •  grad  p) 

—  V  ^  \i  +  grad  u  X  curl  V  —  A  grad  p 
—  curl  curl  p  V  (A7) 

In  order  to  apply  the  boundary. layer  assumptions,  consider  the  x  and  y 
components  of  (A7)  separately.  The  x-wise  or  i  component  in  the  orthogonal 
coordinate  system  is: 


where  3T  is  a  body  force. 
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Now  make  the  usual  order  of  magnitude  analysis  of  (A8).  Let  u#,  p#  and 
x  he  of  order  unity,  that  is,  ue,  p#,  x  -  0(l).  Then  since 


w 


»e 


Also  since, 


Iv^-ow. 


Similarly,  the  magnitude  of  the  other  components  of  (A3)  may  be  found.  They 
are  summarized  in  the  table  below. 


ORDHl 

QUANTITY 

CM 

1 

to 

o 

S2u 

i7 

0(6_1) 

bu 

&  • 

0(1) 

r,  u,  x. 

du  du 
Pt  Sc  * 

0(8) 

bv 

dv  d2v 

0(82) 

v  =  n/p 

d2u  bv 

5F'* 


If  only  terms  of  order  (8)  or  larger —that  is,  terms  of  order  (8),  (l),  (8  ^), 
(5~2) — are  retained,  then  eq. (a8)  reduces  to  the  following  equation. 


fdu  .  „  da  .  „  dul  v  dp  .  u  dr  du  .  du  du  ,  „  d  u 

» [st +  u St  v ■  '‘-rfw'w* 

The  y  or  J  component  of  equation  (A7)  is 


(A9) 
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If  terms  of  order  unity  are  retained,  (AlO)  reduces  to 


p  ay 


(All) 


For  the  flow  under  consideration  there  are  no  body  forces  and  the  flow  is 
steady.  The  momentum  equations  (A9)  and  (All)  then  become 


P 


du 

u  Sic  +  v 


u  dr  du  d 


(A12) 


1 

P 


0(B) 


(A15) 


The  latter  says  that  the  total  change  of  pressure  through  the  boundary  layer 

2 

normal  to  the  surface  is  of  order  B  .  Therefore  it  will  be  neglected,  and 
the  partial  derivative  of  pressure  with  respect  to  x,  dp/dx,  is  replaced 
by  the  total  derivative,  dp/dx,  in  eq.(A12). 


Energy  Equation 

The  basic  form  of  the  energy  equation  for  an  orthogonal  coordinate  system 
is  given  by  Howarth,  page  55/  Reference  5»  In  bis  notation 

+  +  ^(k^-)  (Al4) 

where 

h  =  enthalpy  *  /  c  dT 

J  P 

$  ■  dissipation  function  =  rate  of  work  done  by  viscous  forces 
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and  by  definition 


d 

EC 


b  (* 


(It*7-7) 

(A15) 

dlv  (k  grad  T) 

(Al6) 

The  dissipation  term  In  vector  coordinates  is  given  by 

2  A  9  J  .  /  Tt  V 1 _ A  T7\  A  1?  maAfl  /  J  J  .•  If  \  1  /  a. A  O  \  2  2 


•^2 


J  =  dlv  grad  V  -  2  eLlv(VXcurl  V)  -  2  V*grad(div  V)  +  (curl  V) (div  V)4 

(A17) 

An  order  of  magnitude  analysis  as  was  made  above  for  the  momentum  equation 
gives ,  If  terms  of  only  the  highest  order  are  retained, 

2 


"dh 

St  +  u 


dli  dll  dp  dp  / 
+  v  ^  =  £*-  +  u  ^  +  ti  ( 


+  |,l^)  +r^( 


Sc  5JJ 

The  last  term  can  be  written  In  terms  of  enthalpy  by  the  relation 


(A18) 


.dT  _  k  dh  _  ,  k  dh  u  dh 

Sy  “  cp  Sy  =  M  4Cp  Sy  =  Pr  Sy 

where  Pr  is  the  Prandtl  number.  Substitution  of  (A19)  in  (Al8)  gives 

\2 


(A19) 


Tdh 

St  +  u 


dh  .  _  dhl 

S£  +  v  Sy| 


+  u 


“(If)  +r  &[r“fe!|]  (A20) 


Now  If  steady  flow  is  assumed  and  if  the  momentum  equation  (A12)  is  used  to 
replace  the  term  u  f  (A20)  may  be  written 


[»(S*  ■£)"(£♦*£)] 


d  f  u  f  dh  du  \  / 

Sy  [  Pr  \Sy  Sy  /  +  \ 


„  du\l 

(“ — fe-' 

!“!] 

(  u  ) 

„  du  1 

c  Pr] 

U5y  J 

(A21) 


Now  introduce  total  enthalpy  H 


H  »h+|u2 
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and  substitution  in  (A2l)  gives 


P 


'  dH 

U5T  + 


(A22) 
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APPENDIX  B 


COMPUTER  PROGRAM  FOR  SOLVING  THE  BOUNDARY-LAYER  EQUATIONS 
(IBM  7090  Computer) 

* 

Input 

1.  The  constants: 


He 

V2href 

C1 

href 

a  n 

C2 

>  (Appendix  B) 

n-MAX 

Re  = 

f" 

w 

0 

K  (Section  6.5) 

P-MAX  (Section  6.9) 

L-MAX  (Section  6.9) 

Q-MAX  (Section  6.9) 

(Appendix  B) 

€g  (Appendix  B) 

2 .  The  Table : 

TABLE  1  A 


f 

Not  all  of  the  notation  used  here  are  defined  In  9ectlon  4.0  PRINCIPAL  NOTATION, 
but  If  not^  they  are  defined  when  applicable  here. 
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1.  Determine  local  value  of  PQ  and  Rn  if  they  are  not  inputs. 

(a)  If  n  »  0,  PQ  and  Rq  are  always  inputs 

(b)  If  n  /  0,  calculate  using  (6.13),  (6.14),  (6.97),  (6.98) 

2.  Calculate  Cm  by  (6.26). 


n  -  0. 

<4  =  0,  i  =  0. 


5«  Form  table  of  pe/p  and  C  versus  tj. 

If  n  ■  0,  calculate  gQ 

(a)  If  gw  is  input,  gQ  =  g^  -  (gy  -  l) 


(b)  If  gj  is  input,  gQ  =  (n  -  \)  +  1 

Calculate  pjp  and  C  using  (6.23),  (6.24),  (6.25)  and 

fo  =  Vn. 

n  /  0,  use  table  of  pe/p  and  C  from  n-1  station 

<4  /*  0,  use  values  of  pg/ p  and  C  determined  from  the 
last  solution  of  the  energy  equation. 

7-  Solve  momentua  equation. 

(a)  If  n  =  0,  use  (6.95) 

If  n  ^  0,  use  (6.57) 

If  p  =  0,  use  (6.58)  for  x-wise  derivatives 

If  p  =  1,  use  (6.59)  for  x-vise  derivatives 

(b)  Find  values  of  cp^'  to  be  used  for  the  first  try. 

If  Q  *  0  and  n  =  0,  then  9”  *  f" 

w  wQ 

If  Q  =  0  and  n  /  0,  then  9”  »  q> * ' 

w  wn-l 

If  4  /  0,  then  cp*  *  ■  9" 

W  wft-l 

(c)  Calculate  momentum  equation,  testing 
(l)  If  <p'  >  K,  stop  calculation  and  store  9^' 

as  a  high  solution 
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(2) 

If 

<p*  <  —  K  and  <pM  <0,  stop  calculation 

and  store  q> *  *  as  a  low  solution 

(3) 

If 

n  -  n*  and 

if 

9,(t##)  <  0,  store  as  low  solution 

if 

q)*^)  >0,  store  as  high  solution 

00 

Denote  solution  as  9'  (\j)j_ 

1  and  reenter  at  ^cT^,  using: 

is  high,  then  cp"  =  cp"  —  C. 

wi+l  i  A 

Set 

i  = 

i  + 

(a) 

If 

(O'  • 

wi 

00 

If 

9'' 

i 

is  low ,  then  cp  "  =  cp '  *  +  CL 

wi+l  wi  1 

Repeat  procedures  and  until  at^east  one  high  and 

one  low  value  of  <pT’'  are  Imown.  Enter  at  (  E  J  . 


0- 

10. 


w. 


Reenter  at  ,  using: 

[?w  *J  High  +  fa']  Low 


9" 

i 


where 


Cc] 


High 


=  lowest  high  solution 


=  highest  low  solution 


11. 


Cm  Glow 

Repeat  procedures  and  ©  until  either: 

(a)  there  are  three  trial  solutions  extending  to 

within  the  bounds  of  —  K  <  9'(ti  )  <  K 

00 

Interpolate  by  means  of  (6.43)  and  (6.44)  for  solution 
that  causes  cp1^)  =  0 

00  ['Pw,]High  ~  (3w']low  <  *1*  If  B0>  determin®  the  maxlmum 
value  of  t)  common  to  the  last  three  trial  solutions. 

Denote  this  ij  as  rj-MAX.  Interpolate  by  means  of  (6.43) 
and  (6.44)  for  solution  that  causes  9,(tiw)  a  o,  treating 
t)-MAX  as  tj^. 
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Compare  q  : 


(a)  If  q  *  0,  enter  at  f  f'j 

(b)  If  0  <  q  <  compare 

(1)  If 

>  ‘2 


enter  at  (  F 


(2)  If 


iMrfcli 


enter  at  __ 

(c)  If  q  -  q-MAX,  enter  at  f  CM 


Set  I*  ®  0  • 


(a)  Determine  table  of  C  and  Pr 

(1)  If  q  -  0,  L  =  0,  n  =  0,  use  gQ  from  step  5  to 
determine  C  and  Pr  by  (6.23),  (6.25),  (6.27a), 
(6.27b),  (6.27c) 

(2)  If  q  =  0,  L  *  0,  n  ^  0,  use  C  and  Pr  from 
n-1  station 

(3)  If  q  ji  0,  L  =  0,  use  C  and  Pr  from  q-1  solution 


(k)  If  L  /  0,  use  C  and  Pr  from  L-l  solution 
Solve  energy  equation. 

(a)  If  n  =  0,  use  (6.96) 

If  n  /  0,  use  (6.1+9) 

(b)  Calculation  depends  on  whether  ^  or  ^  is  input 
If  is  input, 

(l)  Calculate,  first  using  =  g^  and  =  gy  —  1. 

Denote  solution  as  ^  ^  1 

(2)  If  ^(t^)  >  0,  solve  for  *2,  -  8^  -  1  —  C2 

(3)  If  ♦!(%„)  <  0,  solve  for  *2,  using  *w  -  gv  “  1  +  Cg 

(1+)  Determine  correct  values  of  ♦  and  if using  (6.56), 
(6.57a),  (6.57b) 
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(c)  If  gw  is  input, 

(1)  Calculate,  first  using  fA  =  g*  •  Denote  solution  as  \|r. 

1  "l  x 

(2)  If  ^(nj  >  0,  solve  for  *2,  using  =  gj  -  c2 

(3)  if  <  °>  solve  for  *2>  1181116  *w  =  S^J  +  c2 

(4)  Determine  correct  values  of  t  and  ♦  using  (6.56), 
(6.57a),  (6.57b) 

Calculate  values  of  C  and  Pr  and  store  as  and  Pr^,  using 

(6.23), (6.25),  (6.27a),  (6.27b),  (6.27c). 

Test  L. 

(a)  If  L  <  L-MAX,  set  L  =  L  +  1  and  enter  at  J 

(b)  If  L  =  L-MAX,  set  the  tables  "" 

CL  =  °Q 
= 

and  solve  for  (pe/p)^  using  (6.24).  Set  q>^'  «  (9^*^ 

Set  Q  =  Q  +  1.  filter  at  f  C  j  . 


Denote  values  of  (C  <p“)S  9",  9',  9 ,  ♦*,  Pe/p,  C,  and  Pf  as 

n  values,  that  is,  (C  9,,),  =  (C  9")^  etc. 

Determine  values  of  f^,  f^,  f^ ' ,  gQ  and  by  (6.36),  (6. h-5 ) * 

^  #  9 

Calculate;  Y  versus  ij  by  (6.106),  by  (6.107),  Q  by  (6. Ill), 
c-  by  (6.116),  and  St  /c_  by  (6.121). 


“J  \  /  J  ******  Vflp  ±  */  \  *  /  » 

Go  to  next  n-station  by  setting  n  =  n  +  1  and  entering  at 
Stop  calculation  when  n  >  n-MAX. 


0  • 
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Print  Out 


1.  Input  constants  and  table. 

2.  The  table 


n 

B 

*w 

a 

e 

a 

°f 

Pr 

w 

3.  For  every  station,  the  table 


n 

f 

D 

D 

(c  <p")' 

s 

s' 

d 

• 

4.  If  "print”  is  indicated  in  table  of  input,  the  following 
is  printed  out  for  every  trial  solution: 

Station  No.  — 


V 

n 

n 

TT* 

C 

Pr 

D 

9 

9’ 

tpi  i 

(C9")i 

C 

BIB 
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04 


(b) 


PRANDTL  NUMBER  AT  WALL,  Prw 


Figure  lO.-Dependence  of  recovery  factor  on  fluid-property  laws. 

(a)  Recovery  factor  versus  Mach  number. 

(b)  Recovery  factor  versus  Prandtl  number. 
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Figure  12. -Variation  of  surface  distance  x,  body  radius  r0,  ue/a*,  and  Mewith  axial  distance  for 
the  rocket  nozzle  problem.  He  =  27.246  x  106ft2/sec2;  gw=  0.17625;  Pr  =  0.78. 
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X 

(b) 


Figure  13.- Results  of  calculation  for  the  rocket  nozzle. 

(a)  Skin  friction,  (b)  Ratio  of  Stanton  number  to  skin  friction. 
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5.0 


AXIAL  DISTANCE,  X 

Figure  14.-Variation  of  local  velocity,  body  radius  and  surface  distance  with  axial  distance 
for  a  spherically  blunted  cone.  M*  =  3.0  and  9.0. 
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Figure  15.-  Continued,  (b)  Skin  friction  and  Stanton  number  for  gw=  0.400. 


Figure  15.- Concluded,  (c)  Skin  friction  and  wall  enthalpy  for  g^,  -  0.0. 
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3.0 


Figure  16.-  Results  of  calculations  for  the  spherically  blunted  cone  at  Mach  number  =  9.0  and  50,000-ft.  altitude, 
(a)  Skin  friction  and  Stanton  number  for  gw=  0.083124. 
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5.0 


Figure  16.- Concluded,  (c)  Skin  friction  and  wall  enthalpy  for  g^=  0.0. 
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Figure  19.- Transverse  curvature  effect  on  shear  parameter  for  flow  over  a  cone  for  Pr  =  1.0,  y  =  1.4. 
Comparison  of  present  method  with  first-order  effect  of  Probstein  and  Elliott. 

Pr  =  1.0;  C  =1.0;  a  =  10°;  Re  =  104/ft. 


(a) 


Figure  20.-  Effect  of  discontinuities  in  gw  on  the  solution  of  flow  over  a  flat  plate  at  Mach  number  =  3.0. 
Pr  =  0.72;  C  =  1.0;  gw=gw  for  x  £1.0;  g„=gW[s  for  x  >  1.0.  Positive  value  of  g^  means 
that  heat  is  being  transferred  to  surface.  Negative  value  means  that  heat  is  being  transferred 
from  surface  to  the  air.  (a)  gw  =  0.90134;  gw  =  0.35714,  0.90134,  and  1.80268. 

W0  wco 
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Figure  21.-  Enthalpy  profiles  for  two  of  the  flows  shown  in  figure  20. 
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s  transfer  on  the  solution  of  flow  over  an  insulated  flat  plate  (g'w=  0.0). 
Arrows  indicate  direction  of  air  flow  out  of  surface. 
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(b) 


Figure  23-  Effect  of  blowing  cold  air  at  the  front  of  a  hot  flat  plate.  External  conditions  and  mass-transfer 
distributions  are  the  same  as  shown  in  figure  22.  Also  variation  of  P  with  x  is  the  same. 

(a)  Wall  enthalpy  gw  versus  x.  (b)  Heat-transfer  parameter  g'w  versus  x. 
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Figure  24.- Effect  of  blowing  hot  air  at  the  front  of  a  cold  flat  plate.  External  conditions  and  mass-transfer 
distributions  are  the  same  as  shown  in  figure  22.  Also  a  variation  of  f*  with  x  is  the  same. 

(a)  Wall  enthalpy  gw  versus  x.  (b)  Heat-transfer  parameter  g'w  versus 
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(c) 


Figure  26.-  Effect  of  blowing  at  the  nose  of  an  insulated  blunt  body  on  the  solution  downstream. 

Mffi  =  3.0;  Pr  =  0.72;  C  =  1.0;  P  =  1.0  for  x  £  1.0;  P  =  0.0  for  x  >  1.0. 

Arrows  indicate  direction  of  air  flow  out  of  surface,  (a)  Pressure  gradient  parameter 
P  versus  x.  (b)  Mass-transfer  distribution,  fw  versus  x.  (c)  Shear  parameter  f'  versus  x. 
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Figure  27 -  Effect  of  blowing  cold  air  at  the  nose  of  a  hot  blunt  body  on  the  solution  downstream. 

External  flow  and  mass-transfer  distributions  are  the  same  as  for  the  flow  in  figure  26. 
(a)  Wall  enthal|}y  gw  versus  x.  (b)  Shear  parameter  f'  versus  x. 
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Concluded. 
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-  (c)  Heat-transfer  p 
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